From 951bf4e9eabcd74d3cef9053fdcb157bbcc70c7c Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 12 Aug 2026 20:38:58 -0400 Subject: [PATCH 1/4] ALL - IMPROVEMENT - Self-describing HDF5 metadata (long_name, units, dims, scales) Every gpec.h5 dataset outside Input/ and the debug-only GalerkinIntegration/Match/ now answers "what is this, in what units, plotted against what" without opening the source, readable natively by h5py/xarray/HDFView: - New Utilities.HDF5Annotations: annotate!(parent, table) applies per-writer path => (; long_name, units, dims) tables post-write (skipping absent paths); make_scale!/attach_scale! wrap the H5DS dimension-scale API (netCDF-4 coordinate mechanism) with Julia-axis -> C-dim index translation; write_root_attrs! stamps schema_version=2.0, Conventions=GPEC-HDF5-2.0, references, title, date_created. - Table-driven, not per-write-call: writers are untouched except one annotate call at the end of each (src/HDF5Schema.jl for the main writer; tables live next to write_galerkin!, the PE writer, KineticForces/Output.jl, and Tearing/Runner/HDF5Output.jl). ~340 dataset annotations total. - Coordinate datasets (psi grids, rational-surface psi, geometry xs/ys) are marked as HDF5 Dimension Scales and attached to the profiles sharing the axis; a greppable "dims" attribute mirrors the scales in Julia axis order. - Attribute wording/units audited by the fortran-physics-reviewer against the layer/field-reconstruction sources; 27 corrections applied (J-weighted field units T*m^2, Q-normalization time vs resistive-kink time, island half-width vs full width, omega_Hz actually rad/s, etc.). - runtests_h5_schema.jl now enforces the contract: long_name+units on every non-exempt dataset, dims on rank >= 2 arrays, root attrs present, scales attached; docs/development/hdf5-conventions.md gains the metadata contract. Attributes are invisible to the rerun leaf-walk and the regression extractor: no tracked value moves. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0129rSTCmYJDBbcH9khHqYnz --- Project.toml | 1 + docs/development/hdf5-conventions.md | 13 + docs/src/utilities.md | 9 + src/ForceFreeStates/Galerkin/GalerkinSolve.jl | 50 +++- src/GeneralizedPerturbedEquilibrium.jl | 7 + src/HDF5Schema.jl | 239 ++++++++++++++++++ src/KineticForces/Output.jl | 53 ++++ src/PerturbedEquilibrium/Utils.jl | 93 +++++++ src/Tearing/Runner/HDF5Output.jl | 104 +++++++- src/Utilities/HDF5Annotations.jl | 93 +++++++ src/Utilities/Utilities.jl | 6 + test/runtests_h5_schema.jl | 43 ++++ 12 files changed, 709 insertions(+), 2 deletions(-) create mode 100644 src/HDF5Schema.jl create mode 100644 src/Utilities/HDF5Annotations.jl 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/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index f62f8777d..b24610d2a 100644 --- a/docs/development/hdf5-conventions.md +++ b/docs/development/hdf5-conventions.md @@ -40,6 +40,19 @@ 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. +- **HDF5 Dimension Scales** (the netCDF-4 coordinate mechanism): shared coordinate datasets (`psi` grids, rational-surface `psi`, geometry `xs`/`ys`) 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)` table next to it (`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). 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/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl index e7e73fd50..c080bf869 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl @@ -240,7 +240,10 @@ 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 + if result.msing == 0 + annotate_galerkin!(out_h5) + return nothing + end out_h5["$gdp/delta"] = result.delta out_h5["$gdp/pest3_A"] = result.Ap out_h5["$gdp/pest3_B"] = result.Bp @@ -287,5 +290,50 @@ function write_galerkin!(out_h5, result::GalerkinResult) end 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/msing" => (; 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"), + "ForceFreeStates/Solutions/GalerkinIntegration/Solution/q" => (; long_name="safety factor on the Galerkin solution grid", dims=("psi",)), + "ForceFreeStates/Solutions/GalerkinIntegration/Solution/issing" => (; long_name="flag: grid node lies on a rational surface", dims=("psi",)), + "ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi" => (; long_name="Galerkin solution functions ξ (arbitrary amplitude)", dims=("mode", "psi", "solution")), + "ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi_deriv" => + (; long_name="ψ_N derivative of the Galerkin solution functions (arbitrary amplitude)", dims=("mode", "psi", "solution")), + "ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi_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" => + (; long_name="outer-region Δ' matrix (2msing×2msing, side-major [L_s1, R_s1, ...]; RDCON Galerkin)", dims=("surface_side", "surface_side")), + "SingularSurfaces/GalerkinDeltaPrime/pest3_A" => (; long_name="PEST-3 matching block A' (Galerkin outer region)", dims=("surface", "surface")), + "SingularSurfaces/GalerkinDeltaPrime/pest3_B" => (; long_name="PEST-3 matching block B' (Galerkin outer region)", dims=("surface", "surface")), + "SingularSurfaces/GalerkinDeltaPrime/pest3_Gamma" => (; long_name="PEST-3 matching block Γ' (Galerkin outer region)", dims=("surface", "surface")), + "SingularSurfaces/GalerkinDeltaPrime/pest3_Delta" => (; long_name="PEST-3 matching block Δ' (Galerkin outer region)", dims=("surface", "surface")), + "SingularSurfaces/GalerkinDeltaPrime/sing_psi" => (; long_name="normalized poloidal flux ψ_N of each rational surface"), + "SingularSurfaces/GalerkinDeltaPrime/sing_q" => (; long_name="safety factor q = m/n at each rational surface", dims=("surface",)), + "SingularSurfaces/GalerkinDeltaPrime/sing_m" => (; long_name="resonant poloidal mode number m at each rational surface", dims=("surface",)), + "SingularSurfaces/GalerkinDeltaPrime/sing_n" => (; long_name="resonant toroidal mode number n at each rational surface", dims=("surface",)), + "SingularSurfaces/GalerkinDeltaPrime/di" => (; long_name="Mercier D_I at each rational surface", dims=("surface",)), + "SingularSurfaces/GalerkinDeltaPrime/alpha" => (; long_name="Frobenius small-solution exponent α at each rational surface", dims=("surface",)), + "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 to everything write_galerkin! wrote. +function annotate_galerkin!(out_h5) + ann = Utilities.HDF5Annotations + ann.annotate!(out_h5, GALERKIN_H5_ANNOTATIONS) + sol = "ForceFreeStates/Solutions/GalerkinIntegration/Solution" + gdp = "SingularSurfaces/GalerkinDeltaPrime" + ann.make_scale!(out_h5, "$sol/psi", "psi") + ann.attach_scale!(out_h5, "$sol/q", 1, "$sol/psi", "psi") + ann.attach_scale!(out_h5, "$sol/issing", 1, "$sol/psi", "psi") + ann.make_scale!(out_h5, "$gdp/sing_psi", "psi_rational") + for a in ("sing_q", "sing_m", "sing_n", "di", "alpha") + ann.attach_scale!(out_h5, "$gdp/$a", 1, "$gdp/sing_psi", "psi_rational") + end return nothing end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 7d3b4ca0e..656316b75 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -68,6 +68,7 @@ const H5_RAW_FORCING = "Input/RawInputs/ForcingTerms" const H5_RAW_COILS = "Input/RawInputs/Coils" const H5_GIT_VERSION = "Info/git_version" +include("HDF5Schema.jl") include("Rerun.jl") # Import ForceFreeStates types and functions needed for main @@ -726,6 +727,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[H5_GIT_VERSION] = git_version @@ -991,6 +995,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..ace08b783 --- /dev/null +++ b/src/HDF5Schema.jl @@ -0,0 +1,239 @@ +# 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. + +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_index", "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/q1lim" => (; long_name="dq/dψ_N at the integration boundary"), + # --- Equilibrium/ scalars (written per-field when set; superset listed) --- + "Equilibrium/ro" => (; long_name="R-coordinate of the magnetic axis", units="m"), + "Equilibrium/zo" => (; long_name="Z-coordinate of the magnetic axis", units="m"), + "Equilibrium/psio" => (; long_name="total poloidal flux difference |ψ_axis - ψ_boundary|", units="Wb/rad"), + "Equilibrium/rsep" => (; long_name="R-coordinates of the plasma boundary", units="m"), + "Equilibrium/zsep" => (; long_name="Z-coordinates of the plasma boundary", units="m"), + "Equilibrium/rext" => (; long_name="R-coordinates of the plasma edge", units="m"), + "Equilibrium/zext" => (; long_name="Z-coordinates of the plasma edge", units="m"), + "Equilibrium/psi0" => (; long_name="normalized poloidal flux at reference location"), + "Equilibrium/b0" => (; long_name="total magnetic field strength at the axis", units="T"), + "Equilibrium/q0" => (; long_name="safety factor at the magnetic axis"), + "Equilibrium/qmin" => (; long_name="minimum safety factor in the plasma"), + "Equilibrium/qmax" => (; long_name="maximum safety factor in the plasma"), + "Equilibrium/qa" => (; long_name="safety factor at the plasma edge"), + "Equilibrium/q95" => (; long_name="safety factor at the 95% flux surface"), + "Equilibrium/qextrema_psi" => (; long_name="normalized poloidal flux at q-profile extrema"), + "Equilibrium/qextrema_q" => (; long_name="safety factor at q-profile extrema"), + "Equilibrium/mextrema" => (; long_name="number of extrema in the q-profile"), + "Equilibrium/rmean" => (; long_name="mean major radius of the plasma", units="m"), + "Equilibrium/amean" => (; long_name="mean minor radius of the plasma", units="m"), + "Equilibrium/aratio" => (; long_name="aspect ratio R0/a"), + "Equilibrium/kappa" => (; long_name="plasma elongation"), + "Equilibrium/delta1" => (; long_name="upper triangularity"), + "Equilibrium/delta2" => (; long_name="lower triangularity"), + "Equilibrium/bt0" => (; long_name="toroidal field at the axis", units="T"), + "Equilibrium/crnt" => (; long_name="plasma current", units="A"), + "Equilibrium/bwall" => (; long_name="toroidal field at the wall", units="T"), + "Equilibrium/betat" => (; long_name="toroidal beta"), + "Equilibrium/betan" => (; long_name="normalized beta β_N"), + "Equilibrium/betap1" => (; long_name="poloidal beta (definition 1)"), + "Equilibrium/betap2" => (; long_name="poloidal beta (definition 2)"), + "Equilibrium/betap3" => (; long_name="poloidal beta (definition 3)"), + "Equilibrium/betaj" => (; long_name="current-weighted beta"), + "Equilibrium/li1" => (; long_name="internal inductance (definition 1)"), + "Equilibrium/li2" => (; long_name="internal inductance (definition 2)"), + "Equilibrium/li3" => (; long_name="internal inductance (definition 3)"), + "Equilibrium/volume" => (; long_name="plasma volume", units="m^3"), + "Equilibrium/bt_sign" => (; long_name="sign of the toroidal field"), + "Equilibrium/psi_norm" => (; long_name="normalized poloidal flux at the axis"), + "Equilibrium/b_norm" => (; long_name="normalized total field strength at the axis"), + "Equilibrium/psi_axis" => (; long_name="poloidal flux at the magnetic axis", units="Wb/rad"), + "Equilibrium/psi_boundary" => (; long_name="poloidal flux at the plasma boundary", units="Wb/rad"), + "Equilibrium/psi_axis_norm" => (; long_name="normalized poloidal flux at the axis"), + "Equilibrium/psi_boundary_norm" => (; long_name="normalized poloidal flux at the boundary"), + "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/verbose" => (; long_name="flag: equilibrium setup ran with verbose output (diagnostic echo)"), + "Equilibrium/diagnose_src" => (; long_name="flag: equilibrium source-data diagnostics were enabled (diagnostic echo)"), + "Equilibrium/diagnose_maxima" => (; long_name="flag: equilibrium extrema diagnostics were enabled (diagnostic echo)"), + # --- Equilibrium/Profiles/ (1-D profiles on the ψ_N grid xs) --- + "Equilibrium/Profiles/xs" => (; long_name="normalized poloidal flux ψ_N profile grid"), + "Equilibrium/Profiles/2piF" => (; long_name="2π F with F = R B_φ the toroidal field function", units="T*m", dims=("psi",)), + "Equilibrium/Profiles/mu0p" => (; long_name="μ0 × plasma pressure", units="T^2", dims=("psi",)), + "Equilibrium/Profiles/dVdpsi" => (; long_name="flux-surface volume derivative dV/dψ_N", units="m^3", dims=("psi",)), + "Equilibrium/Profiles/q" => (; long_name="safety factor profile", dims=("psi",)), + # --- Equilibrium/Geometry/ (2-D flux-coordinate maps on (xs, ys) = (ψ_N, θ/2π)) --- + "Equilibrium/Geometry/xs" => (; long_name="normalized poloidal flux ψ_N geometry grid"), + "Equilibrium/Geometry/ys" => (; long_name="normalized poloidal angle θ/2π geometry grid"), + "Equilibrium/Geometry/rcoords" => (; long_name="squared minor-radius coordinate r² of the working coordinate map", units="m^2", dims=("psi", "theta")), + "Equilibrium/Geometry/offset" => (; long_name="poloidal-angle offset of the working coordinate map (fraction of 2π)", dims=("psi", "theta")), + "Equilibrium/Geometry/nu" => (; long_name="toroidal-angle offset ν = φ − 2πζ of the working coordinate map", units="rad", dims=("psi", "theta")), + "Equilibrium/Geometry/jac" => (; long_name="Jacobian of the (ψ_N, θ, ζ) working coordinates", units="m^3", dims=("psi", "theta")), + # --- LocalStability/ --- + "LocalStability/di" => (; long_name="Mercier ideal interchange criterion D_I", dims=("psi",)), + "LocalStability/dr" => (; long_name="Glasser-Greene-Johnson resistive interchange criterion D_R", dims=("psi",)), + "LocalStability/ballooning_Delta_prime" => (; long_name="high-n ballooning Δ' (distinct from the tearing Δ')", dims=("psi",)), + "LocalStability/psi" => (; long_name="normalized poloidal flux ψ_N of the ballooning α boundary scan"), + "LocalStability/alpha" => (; long_name="experimental normalized pressure gradient α", dims=("psi",)), + "LocalStability/alpha_critical" => (; long_name="critical normalized pressure gradient α for first ballooning stability", dims=("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"), + "ForceFreeStates/Solutions/ForwardIntegration/q" => (; long_name="safety factor at saved solution snapshots", dims=("psi",)), + "ForceFreeStates/Solutions/ForwardIntegration/xi_psi" => (; long_name="fundamental-matrix solutions ξ^ψ (arbitrary amplitude)", dims=("mode", "solution", "psi")), + "ForceFreeStates/Solutions/ForwardIntegration/u2" => + (; long_name="conjugate momenta of the fundamental-matrix solutions (arbitrary amplitude)", dims=("mode", "solution", "psi")), + "ForceFreeStates/Solutions/ForwardIntegration/dxi_psi" => + (; 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")), + "ForceFreeStates/Solutions/ForwardIntegration/crit" => (; long_name="DCON zero-crossing criterion at saved snapshots", dims=("psi",)), + # --- ForceFreeStates/EdgeScan/ (power-normalized (W, N) pencil energies) --- + "ForceFreeStates/EdgeScan/psi" => (; long_name="normalized poloidal flux ψ_N of the edge truncation scan"), + "ForceFreeStates/EdgeScan/q" => (; long_name="safety factor at scan points", dims=("psi",)), + "ForceFreeStates/EdgeScan/total_energy" => (; long_name="power-normalized total energy of the least-stable free-boundary mode (per unit ⟨|ξ|²⟩)", dims=("psi",)), + "ForceFreeStates/EdgeScan/plasma_energy" => (; long_name="power-normalized plasma energy of the least-stable mode (per unit ⟨|ξ|²⟩)", dims=("psi",)), + "ForceFreeStates/EdgeScan/vacuum_energy" => (; long_name="power-normalized vacuum energy of the least-stable mode (per unit ⟨|ξ|²⟩)", dims=("psi",)), + "ForceFreeStates/EdgeScan/vacuum_eigenvalue" => (; long_name="least vacuum eigenvalue of the (W, N) pencil at scan points", dims=("psi",)), + # --- SingularSurfaces/ --- + "SingularSurfaces/msing" => (; long_name="number of rational (singular) surfaces in the domain"), + "SingularSurfaces/psi" => (; long_name="normalized poloidal flux ψ_N of each rational surface"), + "SingularSurfaces/q" => (; long_name="safety factor q = m/n at each rational surface", dims=("surface",)), + "SingularSurfaces/q1" => (; long_name="dq/dψ_N at each rational surface", dims=("surface",)), + "SingularSurfaces/m" => (; long_name="resonant poloidal mode numbers per surface (0-padded)", dims=("surface", "mode")), + "SingularSurfaces/n" => (; long_name="resonant toroidal mode numbers per surface (0-padded)", dims=("surface", "mode")), + "SingularSurfaces/di0" => (; long_name="Mercier D_I evaluated at each rational surface", dims=("surface",)), + "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",)), + "SingularSurfaces/F" => (; long_name="Glasser-Greene-Johnson coefficient F per surface", dims=("surface",)), + "SingularSurfaces/G" => (; long_name="Glasser-Greene-Johnson coefficient G per surface", dims=("surface",)), + "SingularSurfaces/H" => (; long_name="Glasser-Greene-Johnson coefficient H per surface", dims=("surface",)), + "SingularSurfaces/K" => (; long_name="Glasser-Greene-Johnson coefficient K per surface", dims=("surface",)), + "SingularSurfaces/M" => (; long_name="Glasser-Greene-Johnson coefficient M per surface", dims=("surface",)), + "SingularSurfaces/avg_bsq_over_dpsisq" => (; long_name="flux-surface average ⟨B²/|∇ψ_N|²⟩ per surface", units="T^2*m^2", dims=("surface",)), + "SingularSurfaces/avg_bsq" => (; long_name="flux-surface average ⟨B²⟩ per surface", units="T^2", dims=("surface",)), + "SingularSurfaces/p_local" => (; long_name="μ0 × local pressure at each surface", units="T^2", dims=("surface",)), + "SingularSurfaces/p1_local" => (; long_name="μ0 × dp/dψ_N at each surface", units="T^2", dims=("surface",)), + "SingularSurfaces/v1_local" => (; long_name="dV/dψ_N at each surface", units="m^3", dims=("surface",)), + "SingularSurfaces/delta_prime_matrix" => (; long_name="inter-surface Δ' matrix (PEST3 convention, STRIDE BVP with vacuum coupling)", dims=("surface", "surface")), + "SingularSurfaces/delta_prime_raw" => (; long_name="raw 2msing×2msing outer-region D' matrix, side-major ordering [L_s1, R_s1, ...]", dims=("surface_side", "surface_side")), + "SingularSurfaces/delta_coil" => (; long_name="edge coil-response matrix (edge mode × surface-side)", dims=("mode", "surface_side")), + # --- SingularSurfaces/Kinetic/ --- + "SingularSurfaces/Kinetic/kmsing" => (; long_name="number of kinetic singular surfaces (det(F̄) near-zeros)"), + "SingularSurfaces/Kinetic/psi" => (; long_name="normalized poloidal flux ψ_N of kinetic singular surfaces"), + "SingularSurfaces/Kinetic/q" => (; long_name="safety factor at kinetic singular surfaces"), + "SingularSurfaces/Kinetic/q1" => (; long_name="dq/dψ_N at kinetic singular surfaces"), + "SingularSurfaces/Kinetic/scan_psi" => (; long_name="ψ_N grid of the cond(F̄) scan"), + "SingularSurfaces/Kinetic/scan_cond" => (; long_name="condition number of F̄ along the scan"), + "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", "mode")), + "ForceFreeStates/FreeBoundaryStability/W_plasma" => (; long_name="power-normalized plasma energy matrix (per unit ⟨|ξ|²⟩)", dims=("mode", "mode")), + "ForceFreeStates/FreeBoundaryStability/W_vacuum" => (; long_name="power-normalized vacuum energy matrix (per unit ⟨|ξ|²⟩)", dims=("mode", "mode")), + "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"), +] +const _ELM_KINETIC_LETTERS = vcat(_ELM_IDEAL_LETTERS, [("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")], + ["ForceFreeStates/EulerLagrangeMatrices/Ideal/$l" => (; long_name="ideal " * d, dims=("psi", "mode", "mode")) for (l, d) in _ELM_IDEAL_LETTERS], + ["ForceFreeStates/EulerLagrangeMatrices/Kinetic/$l" => (; long_name="kinetic-modified " * d, dims=("psi", "mode", "mode")) 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) + ann = Utilities.HDF5Annotations + ann.annotate!(out_h5, MAIN_H5_ANNOTATIONS) + ann.annotate!(out_h5, ELM_H5_ANNOTATIONS) + + # Coordinate datasets → dimension scales, attached to the profiles sharing the axis. + fwd = "ForceFreeStates/Solutions/ForwardIntegration" + ann.make_scale!(out_h5, "$fwd/psi", "psi") + ann.attach_scale!(out_h5, "$fwd/q", 1, "$fwd/psi", "psi") + ann.attach_scale!(out_h5, "$fwd/crit", 1, "$fwd/psi", "psi") + for a in ("xi_psi", "u2", "dxi_psi") + ann.attach_scale!(out_h5, "$fwd/$a", 3, "$fwd/psi", "psi") + end + + ann.make_scale!(out_h5, "Equilibrium/Profiles/xs", "psi") + for a in ("2piF", "mu0p", "dVdpsi", "q") + ann.attach_scale!(out_h5, "Equilibrium/Profiles/$a", 1, "Equilibrium/Profiles/xs", "psi") + end + + ann.make_scale!(out_h5, "Equilibrium/Geometry/xs", "psi") + ann.make_scale!(out_h5, "Equilibrium/Geometry/ys", "theta") + for a in ("rcoords", "offset", "nu", "jac") + ann.attach_scale!(out_h5, "Equilibrium/Geometry/$a", 1, "Equilibrium/Geometry/xs", "psi") + ann.attach_scale!(out_h5, "Equilibrium/Geometry/$a", 2, "Equilibrium/Geometry/ys", "theta") + end + + ann.make_scale!(out_h5, "SingularSurfaces/psi", "psi_rational") + for a in ("q", "q1", "di0", "E", "F", "G", "H", "K", "M", + "avg_bsq_over_dpsisq", "avg_bsq", "p_local", "p1_local", "v1_local") + ann.attach_scale!(out_h5, "SingularSurfaces/$a", 1, "SingularSurfaces/psi", "psi_rational") + end + + ann.make_scale!(out_h5, "ForceFreeStates/EdgeScan/psi", "psi") + for a in ("q", "total_energy", "plasma_energy", "vacuum_energy", "vacuum_eigenvalue") + ann.attach_scale!(out_h5, "ForceFreeStates/EdgeScan/$a", 1, "ForceFreeStates/EdgeScan/psi", "psi") + end + + elm = "ForceFreeStates/EulerLagrangeMatrices" + ann.make_scale!(out_h5, "$elm/psi", "psi") + for grp in ("Ideal", "Kinetic"), (l, _) in _ELM_KINETIC_LETTERS + ann.attach_scale!(out_h5, "$elm/$grp/$l", 1, "$elm/psi", "psi") + end + + return out_h5 +end diff --git a/src/KineticForces/Output.jl b/src/KineticForces/Output.jl index 502be6b77..40c8c89e7 100644 --- a/src/KineticForces/Output.jl +++ b/src/KineticForces/Output.jl @@ -58,6 +58,59 @@ function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState) mat_g["matrix_$k"] = mat[:, :, k] end end + + # Metadata pass: method tokens are data-driven, so annotate each method group. + for method_name in keys(g) + annotate_kinetic_forces!(g[method_name]) + end +end + +# 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. +const KF_METHOD_H5_ANNOTATIONS = [ + "nn" => (; 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 2n·δW_k", 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"), + "dTdpsi_real" => (; long_name="Re dT_φ/dψ_N torque density at quadrature points", units="N*m", dims=("psi",)), + "dTdpsi_imag" => (; long_name="Im dT_φ/dψ_N (2n·dδW_k/dψ_N energy density) at quadrature points", units="J", dims=("psi",)), + "T_real" => (; long_name="cumulative toroidal torque T_φ(ψ_N) (trapezoidal)", units="N*m", dims=("psi",)), + "T_imag" => (; long_name="cumulative 2n·δW_k(ψ_N) (trapezoidal)", units="J", dims=("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 of each record"), + "EnergyIntegrals/torque_real" => (; long_name="Re of the record's torque contribution", units="N*m"), + "EnergyIntegrals/torque_imag" => (; long_name="Im of the record's torque contribution", units="N*m"), + "EnergyIntegrals/kinetic_energy_real" => (; long_name="Re of the record's kinetic energy contribution", units="J"), + "EnergyIntegrals/kinetic_energy_imag" => (; long_name="Im of the record's kinetic energy contribution", 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_real_all" => (; long_name="Re of the energy-space torque integrand along all trajectories (concatenated)"), + "EnergyIntegrals/integrand_imag_all" => (; long_name="Im of the energy-space torque integrand along all trajectories (concatenated)"), + "EnergyIntegrals/integral_real_all" => (; long_name="Re of the cumulative energy-space integral along all trajectories (concatenated)"), + "EnergyIntegrals/integral_imag_all" => (; long_name="Im of the cumulative energy-space integral along all trajectories (concatenated)"), + "KineticMatrices/matrix_1" => (; long_name="drift-kinetic coefficient matrix 1 of 6 (Logan 2015 Eqs 7.30-7.35)", dims=("mode", "mode")), + "KineticMatrices/matrix_2" => (; long_name="drift-kinetic coefficient matrix 2 of 6 (Logan 2015 Eqs 7.30-7.35)", dims=("mode", "mode")), + "KineticMatrices/matrix_3" => (; long_name="drift-kinetic coefficient matrix 3 of 6 (Logan 2015 Eqs 7.30-7.35)", dims=("mode", "mode")), + "KineticMatrices/matrix_4" => (; long_name="drift-kinetic coefficient matrix 4 of 6 (Logan 2015 Eqs 7.30-7.35)", dims=("mode", "mode")), + "KineticMatrices/matrix_5" => (; long_name="drift-kinetic coefficient matrix 5 of 6 (Logan 2015 Eqs 7.30-7.35)", dims=("mode", "mode")), + "KineticMatrices/matrix_6" => (; long_name="drift-kinetic coefficient matrix 6 of 6 (Logan 2015 Eqs 7.30-7.35)", dims=("mode", "mode")) +] + +# Attach long_name/units/dims + the ψ_N quadrature scale to one method group. +function annotate_kinetic_forces!(method_g) + ann = Utilities.HDF5Annotations + ann.annotate!(method_g, KF_METHOD_H5_ANNOTATIONS) + ann.make_scale!(method_g, "psi", "psi") + for a in ("dTdpsi_real", "dTdpsi_imag", "T_real", "T_imag") + ann.attach_scale!(method_g, a, 1, "psi", "psi") + end + return nothing end """ diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl index 418a1ddd7..322c6cf9c 100644 --- a/src/PerturbedEquilibrium/Utils.jl +++ b/src/PerturbedEquilibrium/Utils.jl @@ -224,5 +224,98 @@ 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", "mode")), + "ResponseMatrices/surface_inductance" => (; long_name="surface inductance L̃ in root-area-weighted field space", dims=("mode", "mode")), + "ResponseMatrices/permeability" => (; long_name="permeability P̃ = Λ̃·L̃⁻¹ in root-area-weighted field space", dims=("mode", "mode")), + "ResponseMatrices/reluctance" => (; long_name="reluctance ϱ̃ in root-area-weighted field space", dims=("mode", "mode")), + "ResponseMatrices/rootarea_to_area_weight_operator" => (; long_name="operator S = Σ/√A at ψ_lim; b̄ = S·b̃", dims=("mode", "mode")), + "ResponseMatrices/surface_area" => (; long_name="control-surface scalar area A = ∮J|∇ψ|dθ; Φ = A·b̄", units="m^2"), + "Response/psi_n" => (; long_name="normalized poloidal flux ψ_N grid shared by the response profiles"), + "Response/xi_psi" => (; long_name="contravariant radial displacement ξ^ψ = ξ·∇ψ_N", dims=("psi", "mode")), + "Response/xi_psi_J" => (; long_name="Jacobian-weighted contravariant radial displacement J·ξ^ψ", units="m^3", dims=("psi", "mode")), + "Response/xi_theta" => (; long_name="Jacobian-weighted contravariant poloidal displacement J·ξ^θ", units="m^3", dims=("psi", "mode")), + "Response/xi_zeta" => (; long_name="Jacobian-weighted contravariant toroidal displacement J·ξ^ζ", units="m^3", dims=("psi", "mode")), + "Response/xi_theta_reg" => (; long_name="regularized Jacobian-weighted contravariant poloidal displacement J·ξ^θ", units="m^3", dims=("psi", "mode")), + "Response/xi_zeta_reg" => (; long_name="regularized Jacobian-weighted contravariant toroidal displacement J·ξ^ζ", units="m^3", dims=("psi", "mode")), + "Response/xi_cova_psi" => (; long_name="covariant radial displacement ξ_ψ", units="m^2", dims=("psi", "mode")), + "Response/xi_cova_theta" => (; long_name="covariant poloidal displacement ξ_θ", units="m^2", dims=("psi", "mode")), + "Response/xi_cova_zeta" => (; long_name="covariant toroidal displacement ξ_ζ", units="m^2", dims=("psi", "mode")), + "Response/clebsch_psi" => (; long_name="Clebsch displacement component ξ^ψ (PENTRC input, gpout_xclebsch convention)", dims=("psi", "mode")), + "Response/clebsch_psi1" => (; long_name="regularized ψ_N derivative of ξ^ψ (× singfac²/(singfac²+reg_spot²))", dims=("psi", "mode")), + "Response/clebsch_alpha" => (; long_name="Clebsch displacement component ξ^α/χ₁ (PENTRC input, gpout_xclebsch convention)", dims=("psi", "mode")), + "Response/xi_n" => (; long_name="physical normal displacement ξ_n", units="m", dims=("psi", "mode")), + "Response/xi_R" => (; long_name="cylindrical displacement component ξ_R (mode space)", units="m", dims=("psi", "mode")), + "Response/xi_Z" => (; long_name="cylindrical displacement component ξ_Z (mode space)", units="m", dims=("psi", "mode")), + "Response/xi_phi" => (; long_name="cylindrical displacement component ξ_φ (mode space)", units="m", dims=("psi", "mode")), + "Response/b_psi_area_weighted" => (; long_name="area-normalized radial field b^ψ/⟨J|∇ψ|⟩_θ", units="T", dims=("psi", "mode")), + "Response/b_n" => (; long_name="physical normal field b_n", units="T", dims=("psi", "mode")), + "Response/b_theta" => (; long_name="Jacobian-weighted contravariant poloidal field J·b^θ", units="T*m^2", dims=("psi", "mode")), + "Response/b_zeta" => (; long_name="Jacobian-weighted contravariant toroidal field J·b^ζ", units="T*m^2", dims=("psi", "mode")), + "Response/b_theta_reg" => (; long_name="regularized Jacobian-weighted contravariant poloidal field J·b^θ", units="T*m^2", dims=("psi", "mode")), + "Response/b_zeta_reg" => (; long_name="regularized Jacobian-weighted contravariant toroidal field J·b^ζ", units="T*m^2", dims=("psi", "mode")), + "Response/b_cova_psi" => (; long_name="covariant radial field b_ψ", units="T*m", dims=("psi", "mode")), + "Response/b_cova_theta" => (; long_name="covariant poloidal field b_θ", units="T*m", dims=("psi", "mode")), + "Response/b_cova_zeta" => (; long_name="covariant toroidal field b_ζ", units="T*m", dims=("psi", "mode")), + "Response/b_R" => (; long_name="cylindrical field component b_R (mode space)", units="T", dims=("psi", "mode")), + "Response/b_Z" => (; long_name="cylindrical field component b_Z (mode space)", units="T", dims=("psi", "mode")), + "Response/b_phi" => (; long_name="cylindrical field component b_φ (mode space)", units="T", dims=("psi", "mode")), + "SingularCoupling/C_resonant_area_weighted_field" => (; long_name="coupling matrix: applied b̃ → resonant area-weighted field b̄^r = Φ^r/A^r", dims=("surface", "mode")), + "SingularCoupling/C_resonant_current" => (; long_name="coupling matrix: applied b̃ → pitch-resonant current", units="A/T", dims=("surface", "mode")), + "SingularCoupling/C_island_width_sq" => (; long_name="coupling matrix: applied b̃ → squared island half-width", units="1/T", dims=("surface", "mode")), + "SingularCoupling/C_penetrated_area_weighted_field" => (; long_name="coupling matrix: applied b̃ → penetrated area-weighted field", dims=("surface", "mode")), + "SingularCoupling/C_delta_prime" => (; long_name="coupling matrix: applied b̃ → forcing-driven Δ'", units="1/T", dims=("surface", "mode")), + "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",)), + "SingularCoupling/resonant_current" => (; long_name="pitch-resonant current per rational surface", units="A", dims=("surface",)), + "SingularCoupling/island_width_sq" => (; long_name="squared island half-width per rational surface (in ψ_N²)", dims=("surface",)), + "SingularCoupling/penetrated_area_weighted_field" => (; long_name="penetrated area-weighted field per rational surface", units="T", dims=("surface",)), + "SingularCoupling/delta_prime" => (; long_name="forcing-driven tearing Δ' per rational surface (Riccati; response to applied forcing)", dims=("surface",)), + "SingularCoupling/forcing_solution_weights" => (; long_name="weights of the forcing solutions in the singular-coupling decomposition", dims=("surface",)), + "SingularCoupling/rational_area" => (; long_name="scalar surface area A^r of each rational surface", units="m^2", dims=("surface",)), + "SingularCoupling/island_half_width" => (; long_name="island half-width per rational surface (in ψ_N)", dims=("surface",)), + "SingularCoupling/chirikov_parameter" => (; long_name="Chirikov overlap parameter: island half-width / half-distance to the neighbouring rational surface", dims=("surface",)), + "SingularCoupling/rational_psi" => (; long_name="normalized poloidal flux ψ_N of each rational surface"), + "SingularCoupling/rational_q" => (; long_name="safety factor q = m/n at each rational surface", dims=("surface",)), + "SingularCoupling/rational_m_res" => (; long_name="resonant poloidal mode number m at each rational surface", dims=("surface",)), + "SingularCoupling/rational_n" => (; long_name="resonant toroidal mode number n at each rational surface", dims=("surface",)), + "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 to the PerturbedEquilibrium group. +function annotate_pe!(pe_group) + ann = Utilities.HDF5Annotations + ann.annotate!(pe_group, PE_H5_ANNOTATIONS) + ann.make_scale!(pe_group, "Response/psi_n", "psi") + for (path, _) in PE_H5_ANNOTATIONS + startswith(path, "Response/") && path != "Response/psi_n" || continue + ann.attach_scale!(pe_group, path, 1, "Response/psi_n", "psi") + end + ann.make_scale!(pe_group, "SingularCoupling/rational_psi", "psi_rational") + for (path, _) in PE_H5_ANNOTATIONS + startswith(path, "SingularCoupling/") && path != "SingularCoupling/rational_psi" || continue + ann.attach_scale!(pe_group, path, 1, "SingularCoupling/rational_psi", "psi_rational") end + return nothing end diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index 2c1ee6611..f35a2c184 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -35,7 +35,10 @@ function write_slayer_hdf5!(parent::Union{HDF5.File,HDF5.Group}, g = create_group(parent, "Tearing") g["enabled"] = Int(result.enabled) - 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,9 +47,108 @@ 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 +# 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/ising" => (; 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/tauk" => (; long_name="Q-normalization time S^(1/3)·τ_H per surface (Q = −τ_k·ω)", units="s", dims=("surface",)), + "PerSurface/tau_r" => (; long_name="resistive diffusion time τ_R per surface", units="s", dims=("surface",)), + "PerSurface/delta_n" => (; 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/dr_val" => (; long_name="resistive interchange D_R = E + F + H² for the critical-Δ formula (auto-derived from GGJ coefficients unless overridden)", dims=("surface",)), + "PerSurface/dgeo_val" => (; 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/dc_tmp" => (; long_name="critical-Δ offset from χ_∥/χ_⊥ matching (Connor-Hastie-Helander 2015 Eq. 59)", dims=("surface",)), + "PerSurface/dc_type" => (; long_name="per-surface D_c prescription label"), + "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/taua" => (; long_name="Alfvén time τ_A per surface", units="s", dims=("surface",)), + "PerSurface/taur" => (; long_name="resistive diffusion time τ_R per surface", units="s", dims=("surface",)), + "PerSurface/v1" => (; long_name="dV/dψ_N at each surface", units="m^3", dims=("surface",)), + "PerSurface/DpMatrix/real" => (; long_name="Re of the full Δ' matrix coupling the rational surfaces", dims=("surface", "surface")), + "PerSurface/DpMatrix/imag" => (; long_name="Im of the full Δ' matrix coupling the rational surfaces", dims=("surface", "surface")), + "Roots/Q_root_real" => (; long_name="Re of the dispersion-root normalized frequency Q (NaN = no root)"), + "Roots/Q_root_imag" => (; long_name="Im of the dispersion-root normalized frequency Q (NaN = no root)"), + "Roots/omega_Hz" => (; long_name="mode rotation angular frequency of each root", units="rad/s"), + "Roots/gamma_Hz" => (; long_name="growth rate of each root", units="1/s"), + "Roots/no_root" => (; long_name="flag: no usable dispersion root found (Q_root is NaN, ω/γ are placeholders)"), + "LayerWidths/ising" => (; 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/dels_db_real" => (; long_name="Re of the dimensionless layer thickness δ_s/d_β", dims=("surface",)), + "LayerWidths/dels_db_imag" => (; long_name="Im of the dimensionless layer thickness δ_s/d_β", dims=("surface",)), + "LayerWidths/delta_s_real" => (; long_name="Re of the complex resistive layer thickness δ_s (Riccati)", dims=("surface",)), + "LayerWidths/delta_s_imag" => (; long_name="Im of the complex resistive layer thickness δ_s (Riccati)", dims=("surface",)), + "LayerWidths/delta_s_m" => (; 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_real" => (; long_name="Re of the concatenated complex entries (rows delimited by offsets)"), + "flat_imag" => (; long_name="Im of the 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_real" => (; long_name="Re of the sampled normalized frequency Q"), + "Q_imag" => (; long_name="Im of the sampled normalized frequency Q"), + "Delta_real" => (; long_name="Re of the inner-layer matching Δ(Q)"), + "Delta_imag" => (; long_name="Im of the 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_real") && ndims(sg["Q_real"]) == 2 + for a in ("Q_real", "Q_imag", "Delta_real", "Delta_imag") + 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}) diff --git a/src/Utilities/HDF5Annotations.jl b/src/Utilities/HDF5Annotations.jl new file mode 100644 index 000000000..416152c8f --- /dev/null +++ b/src/Utilities/HDF5Annotations.jl @@ -0,0 +1,93 @@ +""" + 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` — a tuple of axis names in Julia (column-major) order, axis 1 first, stored +as the greppable string attribute `dims = "(psi, m)"`. Missing paths are skipped. +""" +function annotate!(parent::Union{HDF5.File,HDF5.Group}, table) + 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, ", ") * ")") + end + return parent +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/runtests_h5_schema.jl b/test/runtests_h5_schema.jl index 4ad698b13..3e11f978a 100644 --- a/test/runtests_h5_schema.jl +++ b/test/runtests_h5_schema.jl @@ -32,6 +32,31 @@ function _collect_bad_groups(h5) return bad end +# Metadata contract (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. +_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 + @testset "gpec.h5 schema naming" begin template_dir = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") @@ -61,6 +86,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 From 0ece9c4fe016b262797c2f2684ec9fc5767c37b7 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 12 Aug 2026 21:22:11 -0400 Subject: [PATCH 2/4] ALL - CLEANUP - Wrap annotation-table lines to the 180-char margin Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0129rSTCmYJDBbcH9khHqYnz --- src/ForceFreeStates/Galerkin/GalerkinSolve.jl | 3 ++- src/Tearing/Runner/HDF5Output.jl | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl index c080bf869..8fe08648e 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl @@ -306,7 +306,8 @@ const GALERKIN_H5_ANNOTATIONS = [ (; long_name="ψ_N derivative of the Galerkin solution functions (arbitrary amplitude)", dims=("mode", "psi", "solution")), "ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi_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")), + "ForceFreeStates/Solutions/GalerkinIntegration/Solution/cut_range" => + (; long_name="ψ_N bounds of the excised resonant + extension cells per surface", dims=("surface", "bound")), "SingularSurfaces/GalerkinDeltaPrime/delta" => (; long_name="outer-region Δ' matrix (2msing×2msing, side-major [L_s1, R_s1, ...]; RDCON Galerkin)", dims=("surface_side", "surface_side")), "SingularSurfaces/GalerkinDeltaPrime/pest3_A" => (; long_name="PEST-3 matching block A' (Galerkin outer region)", dims=("surface", "surface")), diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index f35a2c184..efff4ef0d 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -74,7 +74,8 @@ const TEARING_H5_ANNOTATIONS = [ "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/dr_val" => (; long_name="resistive interchange D_R = E + F + H² for the critical-Δ formula (auto-derived from GGJ coefficients unless overridden)", dims=("surface",)), + "PerSurface/dr_val" => + (; long_name="resistive interchange D_R = E + F + H² for the critical-Δ formula (auto-derived from GGJ coefficients unless overridden)", dims=("surface",)), "PerSurface/dgeo_val" => (; 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",)), From 7825f5d7f77d5d7ac0d412cf74d719fbbdd1a316 Mon Sep 17 00:00:00 2001 From: priyanshlunia <40486607+priyanshlunia@users.noreply.github.com> Date: Fri, 14 Aug 2026 16:43:10 -0400 Subject: [PATCH 3/4] ALL - REFACTOR - Rename TOML config variables to descriptive spellings (issue #287) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Rename initials-style and misnomer TOML input keys across all config sections, with warn-and-remap back-compat: old spellings load with a deprecation warning and produce identical control structs until removal after v2.0.0. The use_parallel/use_riccati boolean pair becomes a single `integrator` enum whose remap mirrors the old dispatch order (parallel wins over riccati, both false means serial, old default use_parallel=true). Key renames: | Section | Old | New | |----------------------|----------------------------|--------------------------------------------| | [ForceFreeStates] | use_parallel + use_riccati | integrator = "stride"|"riccati"|"serial" | | [ForceFreeStates] | parallel_threads | integrator_threads | | [ForceFreeStates] | psiedge | dW_edge_scan_start | | [ForceFreeStates] | nstep, diagnose_ca | (deleted — never implemented) | | [Equilibrium] | newq0 (Int) | q0_override (Float64) | | [Equilibrium] | use_galgrid | use_galerkin_grid | | [PerturbedEquilibrium] | reg_spot | regularization_width | | [SLAYER] | dc_type | delta_crit_type | | [SLAYER] | dr_val | delta_crit_D_R | | [SLAYER] | dgeo_val | delta_crit_geo_factor | Value renames (warn-alias for old spellings): | Key | Old | New | |----------------------|------------|-------------------| | [Equilibrium] grid_type | "ldp" | "rational_packed" | | [KineticForces] f0type | "jkp" | "park" | | [SLAYER] delta_crit_type | "rfitzp" | "fitzpatrick" | Mechanics: new _rename_keys!/_rename_value!/_remap_integrator_keys! helpers next to _drop_deprecated_keys!, wired into the Equilibrium/FFS/PE/KF loaders, Rerun.jl (the Input/gpec_toml_raw echo is the only input record post h5 refactor, so old h5 snapshots rerun through the same remap), and slayer_control_from_toml. SLAYER with integrator != "stride" now errors at config time instead of falling back to the stub Δ' at runtime. HDF5 dataset leaf names (Tearing PerSurface/dr_val, dgeo_val, dc_type) deliberately keep the old spellings for schema stability; renaming them is a future schema-bump candidate. Sweep: all example decks, test fixtures, tests, benchmarks, and docs updated to the new spellings (stability.md also fixes its wrong use_parallel=false default claim); pre-commit toml-no-deprecated-keys pattern extended with the old keys; new test/runtests_toml_backcompat.jl covers the remap machinery and old-vs-new ctrl struct equality. Co-Authored-By: Claude Fable 5 --- .pre-commit-config.yaml | 2 +- benchmarks/benchmark_against_fortran_run.jl | 5 +- ...hmark_coil_ForcingTerms_against_fortran.jl | 1 + benchmarks/benchmark_delta_prime_methods.jl | 2 +- .../benchmark_diiid_ideal_ntv_torque.jl | 6 +- .../benchmark_diiid_kinetic_stability.jl | 4 +- benchmarks/benchmark_q_vs_iota_edge.jl | 2 +- benchmarks/benchmark_threads.jl | 19 ++- docs/Project.toml | 3 + docs/development/hdf5-conventions.md | 2 +- docs/src/citations.md | 2 +- docs/src/developer_notes.md | 8 +- docs/src/equilibrium.md | 6 +- docs/src/stability.md | 18 ++- examples/DIIID-like_SLAYER_example/gpec.toml | 10 +- .../gpec.toml | 10 +- .../gpec.toml | 12 +- examples/DIIID-like_ideal_example/gpec.toml | 10 +- examples/LAR_beta_scan/gpec.toml | 6 +- examples/LAR_epsilon_scan/gpec.toml | 6 +- examples/LAR_ideal_match_test/gpec.toml | 6 +- examples/LAR_resistive_match_test/gpec.toml | 6 +- examples/Solovev_ideal_example/gpec.toml | 12 +- examples/Solovev_ideal_example_3D/gpec.toml | 10 +- .../Solovev_ideal_example_multi_n/gpec.toml | 10 +- .../single_n_1/gpec.toml | 6 +- .../single_n_2/gpec.toml | 6 +- .../Solovev_kinetic_NTV_example/gpec.toml | 12 +- .../gpec.toml | 6 +- examples/a10_kinetic_example/gpec.toml | 8 +- src/Analysis/ForceFreeStates.jl | 6 +- src/Equilibrium/DirectEquilibrium.jl | 18 +-- src/Equilibrium/EquilibriumTypes.jl | 26 ++-- src/Equilibrium/InverseEquilibrium.jl | 16 +-- src/ForceFreeStates/EulerLagrange.jl | 38 ++--- src/ForceFreeStates/ForceFreeStatesStructs.jl | 34 ++--- src/ForceFreeStates/Free.jl | 6 +- src/ForceFreeStates/Riccati.jl | 20 +-- src/ForcingTerms/CoilFourier.jl | 4 +- src/GeneralizedPerturbedEquilibrium.jl | 59 +++++++- src/InnerLayer/SLAYER/LayerInputs.jl | 50 +++---- src/InnerLayer/SLAYER/LayerParameters.jl | 58 ++++---- src/KineticForces/Compute.jl | 2 +- src/KineticForces/EnergyIntegration.jl | 10 +- src/KineticForces/KineticForcesStructs.jl | 2 +- .../FieldReconstruction.jl | 32 ++--- .../PerturbedEquilibrium.jl | 2 +- .../PerturbedEquilibriumStructs.jl | 4 +- src/PerturbedEquilibrium/Utils.jl | 2 +- src/Rerun.jl | 1 + src/Tearing/Dispersion/CoupledFullMatch.jl | 2 +- src/Tearing/Runner/Control.jl | 44 ++++-- src/Tearing/Runner/HDF5Output.jl | 9 +- src/Tearing/Runner/run_slayer.jl | 6 +- test/runtests.jl | 1 + test/runtests_equil.jl | 14 +- test/runtests_eulerlagrange.jl | 2 +- test/runtests_parallel_integration.jl | 40 +++--- test/runtests_rerun_from_h5.jl | 2 +- test/runtests_slayer_inputs.jl | 28 ++-- test/runtests_slayer_params.jl | 38 ++--- test/runtests_slayer_runner.jl | 10 +- test/runtests_toml_backcompat.jl | 132 ++++++++++++++++++ test/runtests_vacuum.jl | 2 +- .../gpec.toml | 10 +- .../gpec.toml | 10 +- .../gpec.toml | 6 +- .../gpec.toml | 10 +- .../gpec.toml | 10 +- .../gpec.toml | 6 +- 70 files changed, 598 insertions(+), 390 deletions(-) create mode 100644 test/runtests_toml_backcompat.jl diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 54e179385..506dd06ff 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -69,5 +69,5 @@ repos: - id: toml-no-deprecated-keys name: 'TOML conventions: no deprecated config keys' language: pygrep - entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|power_bp|power_b|power_r|power_rc)\s*=' + entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|power_bp|power_b|power_r|power_rc|nstep|diagnose_ca|use_parallel|use_riccati|parallel_threads|psiedge|newq0|use_galgrid|reg_spot|dc_type|dr_val|dgeo_val)\s*=' files: ^(examples/.*\.toml|test/test_data/.*\.toml)$ diff --git a/benchmarks/benchmark_against_fortran_run.jl b/benchmarks/benchmark_against_fortran_run.jl index 74793ea59..fc73667d0 100644 --- a/benchmarks/benchmark_against_fortran_run.jl +++ b/benchmarks/benchmark_against_fortran_run.jl @@ -149,6 +149,7 @@ function parse_fortran_run(dir::String) eq_file = _find_string(eq_text, "eq_filename"; default="") jac_type = _find_string(eq_text, "jac_type"; default="hamada") grid_type = _find_string(eq_text, "grid_type"; default="ldp") + grid_type == "ldp" && (grid_type = "rational_packed") # map the legacy Fortran value to the Julia spelling psilow = _find_scalar(eq_text, "psilow"; default=1e-4) psihigh = _find_scalar(eq_text, "psihigh"; default=0.993) mpsi = _find_int(eq_text, "mpsi"; default=128) @@ -416,7 +417,7 @@ function write_gpec_toml_coil( # Match the Fortran sas_flag truncation: integration stops at q = qhigh # (= outermost rational q + dmlim). psiedge from dcon.in (=1.0 → no edge dW scan). @printf(io, "qhigh = %.4f\n", qhigh) - @printf(io, "psiedge = %.4f\n", p.psiedge) + @printf(io, "dW_edge_scan_start = %.4f\n", p.psiedge) println(io, "nn_low = $(p.nn)") println(io, "nn_high = $(p.nn)") println(io, "delta_mlow = $(p.delta_mlow)") @@ -476,7 +477,7 @@ function write_gpec_toml_file( println(io, "local_stability_flag = true") println(io, "vac_flag = true") @printf(io, "qlow = %.4f\n", p.qlow) - @printf(io, "psiedge = %.4f\n", p.psiedge) + @printf(io, "dW_edge_scan_start = %.4f\n", p.psiedge) println(io, "nn_low = $(p.nn)") println(io, "nn_high = $(p.nn)") println(io, "delta_mlow = $(p.delta_mlow)") diff --git a/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl b/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl index 22e2b1938..7dbd6b299 100644 --- a/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl +++ b/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl @@ -108,6 +108,7 @@ function parse_fortran_run(dir::String)::FortranRunParams psihigh = _find_scalar(eq_text, "psihigh"; default=0.99) mtheta_eq = _find_int(eq_text, "mtheta"; default=256) grid_type = _find_string(eq_text, "grid_type"; default="ldp") + grid_type == "ldp" && (grid_type = "rational_packed") # map the legacy Fortran value to the Julia spelling # Toroidal mode number and m-range expansion nn = _find_int(dcon_text, "nn"; default=1) diff --git a/benchmarks/benchmark_delta_prime_methods.jl b/benchmarks/benchmark_delta_prime_methods.jl index 4bf77f179..980841387 100644 --- a/benchmarks/benchmark_delta_prime_methods.jl +++ b/benchmarks/benchmark_delta_prime_methods.jl @@ -23,7 +23,7 @@ function setup_and_run_solovev() ex = joinpath(@__DIR__, "..", "test", "test_data", "regression_solovev_ideal_example") inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_riccati"] = true + inputs["ForceFreeStates"]["integrator"] = "riccati" intr = FFS.ForceFreeStatesInternal(; dir_path=ex) ctrl = FFS.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl index 2263e78ce..627ec7546 100644 --- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl +++ b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl @@ -118,12 +118,12 @@ function build_benchmark_tomldir(eq_file::String) eq_filename = "$eq_name" eq_type = "efit" jac_type = "hamada" -grid_type = "ldp" +grid_type = "rational_packed" psilow = 1e-4 psihigh = 0.993 mpsi = 128 mtheta = 256 -newq0 = 0 +q0_override = 0.0 etol = 1e-7 [Wall] @@ -134,7 +134,7 @@ local_stability_flag = true vac_flag = true force_termination = false -psiedge = 1.00 # No edge-scan truncation (dmlim mechanism removed in develop) +dW_edge_scan_start = 1.00 # No edge-scan truncation (dmlim mechanism removed in develop) qlow = 1.02 qhigh = 1e3 sing_start = 0 diff --git a/benchmarks/benchmark_diiid_kinetic_stability.jl b/benchmarks/benchmark_diiid_kinetic_stability.jl index 9ae37b6b3..0624d71b5 100644 --- a/benchmarks/benchmark_diiid_kinetic_stability.jl +++ b/benchmarks/benchmark_diiid_kinetic_stability.jl @@ -129,7 +129,7 @@ psilow = 0.01 psihigh = 0.993 mpsi = 128 mtheta = 256 -newq0 = 0 +q0_override = 0.0 etol = 1e-7 [Wall] @@ -140,7 +140,7 @@ local_stability_flag = true vac_flag = true force_termination = true # Skip PE+KF post-processing — we only need FFS eigenvalues -psiedge = 1.0 # No edge-scan truncation (dmlim mechanism removed in develop) +dW_edge_scan_start = 1.0 # No edge-scan truncation (dmlim mechanism removed in develop) qlow = 1.02 qhigh = 1e3 sing_start = 0 diff --git a/benchmarks/benchmark_q_vs_iota_edge.jl b/benchmarks/benchmark_q_vs_iota_edge.jl index e258fb457..eadb4dd9c 100644 --- a/benchmarks/benchmark_q_vs_iota_edge.jl +++ b/benchmarks/benchmark_q_vs_iota_edge.jl @@ -21,7 +21,7 @@ const EXAMPLE_DIR = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_examp # Dense ldp reference equilibrium: treat its q(ψ) as ground truth function reference_q() _, eq_config, additional_input = GPE.build_inputs_from_toml(EXAMPLE_DIR) - eq_config.grid_type = "ldp" + eq_config.grid_type = "rational_packed" eq_config.mpsi = 1024 equil = GPE.Equilibrium.setup_equilibrium(eq_config, additional_input) return equil, eq_config diff --git a/benchmarks/benchmark_threads.jl b/benchmarks/benchmark_threads.jl index 048c64c66..9e067ca12 100644 --- a/benchmarks/benchmark_threads.jl +++ b/benchmarks/benchmark_threads.jl @@ -1,5 +1,5 @@ # Thread-scaling benchmark for the bidirectional parallel FM integration. -# Runs the Solovev (N=8) and DIIID-like (N=26) examples with use_parallel=true +# Runs the Solovev (N=8) and DIIID-like (N=26) examples with integrator="stride" # across 1, 2, 4, 8 threads and compares against the serial Riccati path. # # Usage (from JPEC_main root): @@ -7,11 +7,10 @@ using GeneralizedPerturbedEquilibrium, TOML, Printf, Statistics -function run_ffs(ex; use_parallel, use_riccati=false) +function run_ffs(ex; integrator) inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_parallel"] = use_parallel - inputs["ForceFreeStates"]["use_riccati"] = use_riccati + inputs["ForceFreeStates"]["integrator"] = integrator inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=ex) ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; @@ -36,17 +35,17 @@ function run_ffs(ex; use_parallel, use_riccati=false) return real(vac.et[1]), intr.numpert_total end -function timed_run(ex; use_parallel, use_riccati=false, nwarm=1, nrep=2) +function timed_run(ex; integrator, nwarm=1, nrep=2) # Warmup for _ in 1:nwarm - run_ffs(ex; use_parallel, use_riccati) + run_ffs(ex; integrator) end # Timed runs times = Float64[] local et1, N for _ in 1:nrep t0 = time() - et1, N = run_ffs(ex; use_parallel, use_riccati) + et1, N = run_ffs(ex; integrator) push!(times, time() - t0) end return mean(times), et1, N @@ -60,9 +59,9 @@ diiid_ex = joinpath(root, "examples", "DIIID-like_ideal_example") println("\n=== Thread-scaling benchmark ($(nthreads) thread(s)) ===\n") for (label, ex) in [("Solovev", sol_ex), ("DIIID-like", diiid_ex)] - t_std, et_std, N = timed_run(ex; use_parallel=false, use_riccati=false) - t_ric, et_ric, _ = timed_run(ex; use_parallel=false, use_riccati=true) - t_par, et_par, _ = timed_run(ex; use_parallel=true, use_riccati=false) + t_std, et_std, N = timed_run(ex; integrator="serial") + t_ric, et_ric, _ = timed_run(ex; integrator="riccati") + t_par, et_par, _ = timed_run(ex; integrator="stride") err_ric = abs(et_ric - et_std) / abs(et_std) * 100 err_par = abs(et_par - et_std) / abs(et_std) * 100 diff --git a/docs/Project.toml b/docs/Project.toml index 12f07f702..83d3d84a8 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,5 +2,8 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" GeneralizedPerturbedEquilibrium = "462872dd-e066-4d2e-b993-6468b5239634" +[sources] +GeneralizedPerturbedEquilibrium = {path = ".."} + [compat] Documenter = "1.0" diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index b24610d2a..6eae57642 100644 --- a/docs/development/hdf5-conventions.md +++ b/docs/development/hdf5-conventions.md @@ -38,7 +38,7 @@ Top level (10 groups): | `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}`, `Scan/Surface_/` | | `SurfaceGeometries/` | `{Plasma,Wall}/{x,y,z}` point clouds | -Reserved (documented, not yet written): `ForceFreeStates/Solutions/RiccatiIntegration/` — the third integrator backend slot alongside `ForwardIntegration` and `GalerkinIntegration`. +Reserved (documented, not yet written): `ForceFreeStates/Solutions/RiccatiIntegration/` — the third integrator backend slot alongside `ForwardIntegration` and `GalerkinIntegration`, matching the `[ForceFreeStates]` `integrator = "riccati"` option. ## Metadata contract (self-describing datasets) diff --git a/docs/src/citations.md b/docs/src/citations.md index 143b9a493..ea5ddaee3 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -36,7 +36,7 @@ The primary reference for the `ForceFreeStates` module. Derives the Euler-Lagran > *Physics of Plasmas* **25**, 032507 (2018). > DOI: [10.1063/1.5007042](https://doi.org/10.1063/1.5007042) -Reformulates the DCON eigenvalue problem as a Riccati matrix ODE, enabling parallel integration across singular surfaces and faster computation. Implemented in `src/ForceFreeStates/Riccati.jl` and enabled via `use_riccati = true` in `[ForceFreeStates]`. +Reformulates the DCON eigenvalue problem as a Riccati matrix ODE, enabling parallel integration across singular surfaces and faster computation. Implemented in `src/ForceFreeStates/Riccati.jl` and enabled via `integrator = "riccati"` in `[ForceFreeStates]`. --- diff --git a/docs/src/developer_notes.md b/docs/src/developer_notes.md index 885a3de3c..aad8640f7 100644 --- a/docs/src/developer_notes.md +++ b/docs/src/developer_notes.md @@ -69,7 +69,7 @@ Three things follow. First, convergence is per-surface: the outer surface to about 1% — while `dpm[2,2]` is marginal and `dpm[1,1]` (q=2) never settles, still moving ~7% between the two tightest grids. A plateau detector must therefore report per-surface rather than pass/fail for the whole diagonal. -(As the `ldp` scan below shows, q=2 is not inherently unconvergeable — it is the +(As the `rational_packed` (formerly `ldp`) scan below shows, q=2 is not inherently unconvergeable — it is the auto grid that prevents it from settling.) Second, the growth rate is linear in Δ′: `gamma/dpm[1,1]` is 24.1 to within 0.5% @@ -83,7 +83,7 @@ after pass 2, so tightening `psi_accuracy` moves it *further* from self-consistency rather than closer, and the warning's advice to "consider tightening psi_accuracy" is counterproductive in this regime. -The same deck on a deterministic `ldp` grid, which skips the measure-and-re-form +The same deck on a deterministic `rational_packed` grid, which skips the measure-and-re-form step entirely, converges: | mpsi | dpm[1,1] | dpm[2,2] | dpm[3,3] | gamma 2/1 (Hz) | @@ -102,14 +102,14 @@ value should do. The non-convergence under the auto grid is an artifact of the generator, not of the Δ′ extraction. Two consequences. A plateau criterion is implementable today against a fixed -`ldp` grid, without waiting on the auto-grid work. And the auto grid's answers +`rational_packed` grid, without waiting on the auto-grid work. And the auto grid's answers are biased in both directions relative to the converged value: at its default `psi_accuracy` it gave 6.39 (28% low), at its tightest 9.51 (7% high). Anything pinned on the auto grid should be read with that in mind. This is not implemented. Doing it properly needs: - - either a fixed `ldp` grid (which already converges, see above) or knot + - either a fixed `rational_packed` grid (which already converges, see above) or knot refinement iterated to a fixed point (repeat the measure-and-re-form step until `implied_knot_count` stops exceeding the grid in use), since without it the scan target keeps moving; diff --git a/docs/src/equilibrium.md b/docs/src/equilibrium.md index ecff97671..056e99192 100644 --- a/docs/src/equilibrium.md +++ b/docs/src/equilibrium.md @@ -96,18 +96,18 @@ built by a **two-pass measured-curvature refinement** driven by the single accur (no file re-read). Every region's knot count scales as τ^(-1/3), so tightening `psi_accuracy` refines -the core, pedestal, and edge proportionally. The legacy `grid_type = "ldp"` +the core, pedestal, and edge proportionally. The legacy `grid_type = "rational_packed"` (formerly `"ldp"`) (sin²-spaced), `"pow1"`, `"uniform"`, and explicit `mpsi > 0` (single-pass, fixed layout) are still supported. Library users calling `setup_equilibrium` directly with `mpsi = 0` receive the coarse pass-1 grid; use `refined_psi_grid` and the `override_psi_nodes` keyword to apply the refinement manually. -The packing on the DIII-D-like example (n=1) compared to fixed `ldp` grids — note the +The packing on the DIII-D-like example (n=1) compared to fixed `rational_packed` grids — note the coarse spacing across the smooth mid-radius, the spacing dips at each rational surface, and the core/pedestal/edge packing (`benchmarks/plot_grid_knot_placement.jl` regenerates this figure): -![Radial knot packing: auto two-pass vs ldp](assets/grid_knot_placement.png) +![Radial knot packing: auto two-pass vs rational_packed](assets/grid_knot_placement.png) Decomposing the density by source on the same example shows the pedestal band (ψ_N ≈ 0.85–0.98) is driven by *measured* curvature, not the edge floor: the pressure, diff --git a/docs/src/stability.md b/docs/src/stability.md index 3e0474482..54234a1a6 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -55,12 +55,12 @@ numerical strategies. columns of ``U_2`` that correspond to resonant modes are zeroed via Gaussian reduction (GR), keeping the solution bounded. This is the reference path for correctness comparisons. -Enable with (default): +Enable with: ```toml [ForceFreeStates] -use_riccati = false -use_parallel = false +integrator = "serial" ``` +(The default is `integrator = "stride"`, the fundamental-matrix path below.) ### Riccati integration @@ -87,8 +87,7 @@ directly in column `ipert_res` — without Gaussian reduction — and renormaliz Enable with: ```toml [ForceFreeStates] -use_riccati = true -use_parallel = false +integrator = "riccati" ``` **Speedup** (benchmarked on reference examples): @@ -126,10 +125,10 @@ The implementation uses a `direction` field on `IntegrationChunk`: crossing chunk. `balance_integration_chunks` preserves this: the sub-chunk closest to the rational surface inherits `direction`, while the earlier sub-chunk always gets `direction=+1`. -Enable with: +Enable with (default): ```toml [ForceFreeStates] -use_parallel = true +integrator = "stride" ``` **Accuracy** (N=26, DIIID-like example): energy eigenvalue within 2% of standard path. @@ -240,8 +239,7 @@ All `ForceFreeStates` options are set in the `[ForceFreeStates]` section of `gpe ```toml [ForceFreeStates] # Integration driver -use_riccati = false # true: Riccati path (faster, same accuracy) -use_parallel = false # true: parallel FM path (multi-thread, large N) +integrator = "stride" # "stride" (default): parallel FM path (multi-thread, Δ' matrix); "riccati": Riccati path; "serial": reference shooting path # Mode space nn_low = 1 # lowest toroidal mode number @@ -307,7 +305,7 @@ metric = FFS.make_metric(equil, intr.mpert) ffit = FFS.make_matrix(equil, intr, metric) # Choose integration driver. The top-level `eulerlagrange_integration` dispatches -# to the parallel or Riccati path based on ctrl.use_parallel / ctrl.use_riccati, +# on ctrl.integrator ("stride", "riccati", or "serial"), # and always returns a 4-tuple (odet, propagators, chunks, S_at_surface_left). odet, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr) diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index a89b382f6..300f3bcd1 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -14,7 +14,7 @@ psihigh = 0.9995 # Upper limit of normalized flux coordinate mpsi = 0 # Number of radial grid points (0 = auto-compute from psi_accuracy) psi_accuracy = 0.001 # Target absolute error in q for auto-mpsi mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -33,7 +33,7 @@ force_termination = true # Run FFS + SLAYER, skip PerturbedEquilibrium local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW scan band: diagnostic dW(ψ) computed for ψ ∈ [psiedge, psilim]; integration domain set by qhigh / psihigh / dmlim +dW_edge_scan_start = 0.99 # Edge dW scan band: diagnostic dW(ψ) computed for ψ ∈ [dW_edge_scan_start, psilim]; integration domain set by qhigh / psihigh / dmlim qlow = 1.02 # Integration initiated at q determined by min(q0, qlow)... qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh)... sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) @@ -52,8 +52,8 @@ singfac_min = 1e-4 # Fractional distance from rational q at which id ucrit = 1e4 # Maximum fraction of solutions allowed before re-normalized # Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 1 # serial/bit-deterministic BVP — keeps the regression Δ' (and hence γ) reproducible +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (unlocks SingularSurfaces/delta_prime_matrix) +integrator_threads = 1 # serial/bit-deterministic BVP — keeps the regression Δ' (and hence γ) reproducible populate_dense_xi = false # No PerturbedEquilibrium here; the dense EL pass has no consumer (auto-disabled under force_termination anyway). SLAYER needs only delta_prime_matrix from the parallel BVP. set_psilim_via_dmlim = true # TRUE for diverted geqdsks — q → ∞ at separatrix, so dmlim truncation avoids the δW kink instability at negligible domain cost dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim) / n @@ -69,7 +69,7 @@ enabled = true # Run the SLAYER tearing-mode analysis inner_model = "slayer_fitzpatrick" # Inner-layer Δ(Q) model: "slayer_fitzpatrick", "ggj_shooting", or "ggj_galerkin" scan_mode = "amr" # Q-plane scan strategy: "amr" (adaptive refinement) or "brute_force" coupling_mode = "uncoupled" # "uncoupled" (per-surface) or "coupled" (multi-surface determinant) -dc_type = "none" # Critical-Δ offset selector: "none", "lar", "rfitzp", or "toroidal" +delta_crit_type = "none" # Critical-Δ offset selector: "none", "lar", "fitzpatrick", or "toroidal" mu_i = 2.0 # Ion mass in proton-mass units (2.0 = deuterium) zeff = 1.0 # Effective charge chi_perp = 1.0 # fallback only; the kinetic file supplies χ⊥(ψ) diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml index f84d8c52a..f017965b9 100644 --- a/examples/DIIID-like_gal_resistive_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_example/gpec.toml @@ -8,13 +8,13 @@ eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file eq_type = "efit" # Type of the input 2D equilibrium file jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type (ldp = linear-derivative packing toward rationals) +grid_type = "rational_packed" # Radial grid packing type (rational_packed = sin²-packed radial grid) psilow = 1e-4 # Lower limit of normalized flux coordinate psihigh = 0.993 # Upper limit of normalized flux coordinate (0.993 stays clear of the separatrix; truncating at ≳0.998 is numerically unreasonable here) mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) psi_accuracy = 0.001 # Target relative accuracy of splined profile derivatives for the auto grid mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -32,7 +32,7 @@ equal_arc_wall = true # Equal arc length distribution of nodes on wall local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW scan band: diagnostic dW(ψ) computed for ψ ∈ [psiedge, psilim]; integration domain set by qhigh / psihigh / dmlim +dW_edge_scan_start = 0.99 # Edge dW scan band: diagnostic dW(ψ) computed for ψ ∈ [dW_edge_scan_start, psilim]; integration domain set by qhigh / psihigh / dmlim qlow = 1.02 # Integration initiated at q determined by min(q0, qlow)... qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh)... sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) @@ -51,8 +51,8 @@ singfac_min = 1e-4 # Fractional distance from rational q at which id ucrit = 1e4 # Maximum fraction of solutions allowed before re-normalized # Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -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) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (unlocks SingularSurfaces/delta_prime_matrix) +integrator_threads = 2 # Stride BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) populate_dense_xi = false # Dense axis-basis ξ for the FFS HDF5 output. Not needed here: no [PerturbedEquilibrium] section, and the gal-matched path builds its own dense ξ. Set true only for a shooting-fed PE run. set_psilim_via_dmlim = false # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim) / n (only used when set_psilim_via_dmlim = true) diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml index e5536c75e..e37976074 100644 --- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml @@ -7,13 +7,13 @@ eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file eq_type = "efit" # Type of the input 2D equilibrium file jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type (ldp = linear-derivative packing toward rationals) +grid_type = "rational_packed" # Radial grid packing type (rational_packed = sin²-packed radial grid) psilow = 1e-4 # Lower limit of normalized flux coordinate psihigh = 0.993 # Upper limit of normalized flux coordinate (0.993 stays clear of the separatrix; truncating at ≳0.998 is numerically unreasonable here) mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) psi_accuracy = 0.001 # Target relative accuracy of splined profile derivatives for the auto grid mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -31,7 +31,7 @@ equal_arc_wall = true # Equal arc length distribution of nodes on wall local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW scan band: diagnostic dW(ψ) computed for ψ ∈ [psiedge, psilim]; integration domain set by qhigh / psihigh / dmlim +dW_edge_scan_start = 0.99 # Edge dW scan band: diagnostic dW(ψ) computed for ψ ∈ [dW_edge_scan_start, psilim]; integration domain set by qhigh / psihigh / dmlim qlow = 1.02 # Integration initiated at q determined by min(q0, qlow)... qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh)... sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) @@ -50,8 +50,8 @@ singfac_min = 1e-4 # Fractional distance from rational q at which id ucrit = 1e4 # Maximum fraction of solutions allowed before re-normalized # Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -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) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (unlocks SingularSurfaces/delta_prime_matrix) +integrator_threads = 2 # Stride BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) populate_dense_xi = false # Dense axis-basis ξ for the FFS HDF5 output. Not needed here: no [PerturbedEquilibrium] section, and the gal-matched path builds its own dense ξ. Set true only for a shooting-fed PE run. set_psilim_via_dmlim = false # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim) / n (only used when set_psilim_via_dmlim = true) @@ -100,4 +100,4 @@ compute_response = true # Compute plasma response to forcing compute_singular_coupling = true # Compute singular layer coupling metrics verbose = true # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 -reg_spot = 0.05 # Regularization width for singular surfaces (0 = disabled) +regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 9a3bf9bc5..93a4da7ec 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -13,7 +13,7 @@ psihigh = 0.995 # Upper limit of normalized poloidal flux (captur mpsi = 0 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) psi_accuracy = 0.001 # Target relative accuracy of splined profile derivatives for the auto grid mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -31,7 +31,7 @@ equal_arc_wall = true # Equal arc length distribution of nodes on wall local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) @@ -50,8 +50,8 @@ singfac_min = 1e-4 # Fractional distance from rational q at which id ucrit = 1e4 # 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) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (unlocks SingularSurfaces/delta_prime_matrix) +integrator_threads = 2 # Stride 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 = true # Truncate at (last_rational_q + dmlim)/n — TRUE for diverted equilibria (q → ∞ at separatrix) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) @@ -73,7 +73,7 @@ compute_response = true # Compute plasma response to forcing compute_singular_coupling = true # Compute singular layer coupling metrics verbose = true # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 -reg_spot = 0.05 # Regularization width for singular surfaces (0 = disabled) +regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) [KineticForces] kinetic_file = "TkMkr_D3Dlike_Hmode_kinetic.h5" # GPEC HDF5 kinetic schema (psi,n_i,n_e,T_i,T_e,omega_E,+chi_e,chi_phi). Legacy .gpeckf/.kin still readable. diff --git a/examples/LAR_beta_scan/gpec.toml b/examples/LAR_beta_scan/gpec.toml index 8ac649d06..d75d46e1f 100644 --- a/examples/LAR_beta_scan/gpec.toml +++ b/examples/LAR_beta_scan/gpec.toml @@ -10,7 +10,7 @@ [Equilibrium] eq_type = "tj_analytic" # TJ-analytic model (inverse pipeline; Fitzpatrick https://github.com/rfitzp/TJ) jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 0.01 # Lower limit of normalized poloidal flux psihigh = 0.995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) @@ -52,8 +52,8 @@ singfac_min = 1e-4 # Fractional distance from rational q at which id ucrit = 1e4 # Column-norm threshold that triggers solution renormalization sing_order = 6 # Order of the singular-surface (Frobenius) series expansion -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) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (unlocks SingularSurfaces/delta_prime_matrix) +integrator_threads = 2 # Stride BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) populate_dense_xi = false # Append serial-EL pass for dense ξ; not needed without [PerturbedEquilibrium] (default false) set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/examples/LAR_epsilon_scan/gpec.toml b/examples/LAR_epsilon_scan/gpec.toml index f0d174ecb..e96f3de2b 100644 --- a/examples/LAR_epsilon_scan/gpec.toml +++ b/examples/LAR_epsilon_scan/gpec.toml @@ -10,7 +10,7 @@ [Equilibrium] eq_type = "tj_analytic" # TJ-analytic model (inverse pipeline; overridden to "tj_analytic_direct" by run_scan.jl) jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 0.01 # Lower limit of normalized poloidal flux psihigh = 0.995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) @@ -53,8 +53,8 @@ singfac_min = 1e-4 # Fractional distance from rational q at which id ucrit = 1e4 # Column-norm threshold that triggers solution renormalization sing_order = 6 # Order of the singular-surface (Frobenius) series expansion -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) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (unlocks SingularSurfaces/delta_prime_matrix) +integrator_threads = 2 # Stride BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) populate_dense_xi = false # Append serial-EL pass for dense ξ; not needed without [PerturbedEquilibrium] (default false) set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/examples/LAR_ideal_match_test/gpec.toml b/examples/LAR_ideal_match_test/gpec.toml index d150cbac2..db60498f8 100644 --- a/examples/LAR_ideal_match_test/gpec.toml +++ b/examples/LAR_ideal_match_test/gpec.toml @@ -9,7 +9,7 @@ [Equilibrium] eq_type = "tj_analytic" # Type of the input 2D equilibrium (analytic large-aspect-ratio model) jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type (ldp = linear-derivative packing toward rationals) +grid_type = "rational_packed" # Radial grid packing type (rational_packed = sin²-packed radial grid) psilow = 0.01 # Lower limit of normalized flux coordinate psihigh = 0.995 # Upper limit of normalized flux coordinate mpsi = 128 # Number of radial grid points (0 = auto-compute from psi_accuracy) @@ -50,8 +50,8 @@ ucrit = 1e4 # Maximum fraction of solutions allowed before r sing_order = 6 # Power-series order for the singular-surface asymptotics save_interval = 3 # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. -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 is about +20% speedup) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (unlocks SingularSurfaces/delta_prime_matrix) +integrator_threads = 2 # Stride BVP thread cap (1 = serial/bit-deterministic; 2 is about +20% speedup) populate_dense_xi = false # Dense axis-basis xi for the FFS HDF5 output. Not needed here: the gal-matched path builds its own dense xi. set_psilim_via_dmlim = false # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim) / n (only used when set_psilim_via_dmlim = true) diff --git a/examples/LAR_resistive_match_test/gpec.toml b/examples/LAR_resistive_match_test/gpec.toml index 9900ad8e7..c02191605 100644 --- a/examples/LAR_resistive_match_test/gpec.toml +++ b/examples/LAR_resistive_match_test/gpec.toml @@ -10,7 +10,7 @@ [Equilibrium] eq_type = "tj_analytic" # Type of the input 2D equilibrium (analytic large-aspect-ratio model) jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type (ldp = linear-derivative packing toward rationals) +grid_type = "rational_packed" # Radial grid packing type (rational_packed = sin²-packed radial grid) psilow = 0.01 # Lower limit of normalized flux coordinate psihigh = 0.995 # Upper limit of normalized flux coordinate mpsi = 128 # Number of radial grid points (0 = auto-compute from psi_accuracy) @@ -51,8 +51,8 @@ ucrit = 1e4 # Maximum fraction of solutions allowed before r sing_order = 6 # Power-series order for the singular-surface asymptotics save_interval = 3 # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. -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 is about +20% speedup) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (unlocks SingularSurfaces/delta_prime_matrix) +integrator_threads = 2 # Stride BVP thread cap (1 = serial/bit-deterministic; 2 is about +20% speedup) populate_dense_xi = false # Dense axis-basis xi for the FFS HDF5 output. Not needed here: the gal-matched path builds its own dense xi. set_psilim_via_dmlim = false # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim) / n (only used when set_psilim_via_dmlim = true) diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml index 66ba0d48b..edcfff2c6 100644 --- a/examples/Solovev_ideal_example/gpec.toml +++ b/examples/Solovev_ideal_example/gpec.toml @@ -6,12 +6,12 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -42,7 +42,7 @@ compute_response = true # Compute plasma response to forcing compute_singular_coupling = true # Compute singular layer coupling metrics verbose = true # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 -reg_spot = 0.05 # Regularization width for singular surfaces (0 = disabled) +regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) # Note: SLAYER tearing analysis is not run on the Solovev analytic equilibrium # — its Δ' / inner-layer dispersion does not yield meaningful tearing roots. @@ -53,7 +53,7 @@ reg_spot = 0.05 # Regularization width for singular surf local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) @@ -72,8 +72,8 @@ ucrit = 1e3 # Column-norm threshold that triggers solution ren save_interval = 3 # Save every Nth ODE step (1=all). Always saves near rational surfaces. # Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -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) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (unlocks SingularSurfaces/delta_prime_matrix) +integrator_threads = 2 # Stride 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 dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml index 70123040b..1fcf67efb 100644 --- a/examples/Solovev_ideal_example_3D/gpec.toml +++ b/examples/Solovev_ideal_example_3D/gpec.toml @@ -5,12 +5,12 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -18,7 +18,7 @@ force_termination = false # Terminate after equilibrium setup (ski local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 1.0 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +dW_edge_scan_start = 1.0 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) @@ -38,8 +38,8 @@ ucrit = 1e3 # Column-norm threshold that triggers solution ren save_interval = 3 # Save every Nth ODE step (1=all). Always saves near rational surfaces. # Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -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) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (unlocks SingularSurfaces/delta_prime_matrix) +integrator_threads = 2 # Stride BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) populate_dense_xi = false # Append serial-EL pass for dense ξ; not needed without [PerturbedEquilibrium] (default false) set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/examples/Solovev_ideal_example_multi_n/gpec.toml b/examples/Solovev_ideal_example_multi_n/gpec.toml index 3fb3f7b06..1a52f1a93 100644 --- a/examples/Solovev_ideal_example_multi_n/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/gpec.toml @@ -6,12 +6,12 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -30,7 +30,7 @@ equal_arc_wall = true # Equal arc length distribution of nodes local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) @@ -48,8 +48,8 @@ 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 (multi-n Δ' matrix has open issues — sing_lim! warns and skips — but ξ and energies are valid) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (multi-n Δ' matrix has open issues — sing_lim! warns and skips — but ξ and energies are valid) +integrator_threads = 2 # Stride BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) populate_dense_xi = false # Append serial-EL pass for dense ξ; not needed without [PerturbedEquilibrium] (default false) set_psilim_via_dmlim = false # FALSE for multi-n — dmlim truncation is ambiguous when n varies (sing_lim! skips anyway) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml index 8b2b75007..cc34dfeeb 100644 --- a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml @@ -5,12 +5,12 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -29,7 +29,7 @@ equal_arc_wall = true # Equal arc length distribution of nodes local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) diff --git a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml index c7a51c723..e20a0145d 100644 --- a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml @@ -5,12 +5,12 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -29,7 +29,7 @@ equal_arc_wall = true # Equal arc length distribution of nodes local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index 733a92c31..22d95f19e 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/examples/Solovev_kinetic_NTV_example/gpec.toml @@ -8,12 +8,12 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, other) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid points (0 = auto-compute from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -44,13 +44,13 @@ compute_response = true # Compute plasma response to forcing compute_singular_coupling = true # Compute singular layer coupling metrics verbose = true # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 -reg_spot = 0.05 # Regularization width for singular surfaces (0 = disabled) +regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) [ForceFreeStates] local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) @@ -69,8 +69,8 @@ ucrit = 1e3 # Column-norm threshold that triggers solution ren save_interval = 3 # Save every Nth ODE step (1=all). Always saves near rational surfaces. # Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -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) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (unlocks SingularSurfaces/delta_prime_matrix) +integrator_threads = 2 # Stride 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 dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml index 149f40af7..0b2e1dab0 100644 --- a/examples/Solovev_kinetic_calculated_example/gpec.toml +++ b/examples/Solovev_kinetic_calculated_example/gpec.toml @@ -6,12 +6,12 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized flux coordinate psihigh = 0.9995 # Upper limit of normalized flux coordinate mpsi = 16 # Number of radial grid points mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -29,7 +29,7 @@ equal_arc_wall = true # Equal arc length distribution of nodes local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW scan band: dW(ψ) computed for ψ ∈ [psiedge, psilim], integration truncated at peak +dW_edge_scan_start = 0.99 # Edge dW scan band: dW(ψ) computed for ψ ∈ [dW_edge_scan_start, psilim], integration truncated at peak qlow = 1.02 # Integration initiated at q determined by min(q0, qlow)... qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh)... sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index d3f61b391..a2855eaaf 100644 --- a/examples/a10_kinetic_example/gpec.toml +++ b/examples/a10_kinetic_example/gpec.toml @@ -10,12 +10,12 @@ eq_filename = "fix_a100_k10_q2_bn010_prof1" # Path to equilibrium file eq_type = "efit" # Type of the input 2D equilibrium file jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-3 # Lower limit of normalized poloidal flux psihigh = 0.99 # Upper limit of normalized poloidal flux mpsi = 16 # Number of radial grid points (low value for fast iteration) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver [Wall] @@ -26,7 +26,7 @@ local_stability_flag = true # Perform local stability analysis (Mercier and b vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge -psiedge = 1.0 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +dW_edge_scan_start = 1.0 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) @@ -55,6 +55,6 @@ zimp = 6 # Impurity charge mimp = 12 # Impurity mass electron = false # Include electron contribution (false = ion-only) nutype = "harmonic" # Collision operator (zero, small, krook, harmonic) -f0type = "maxwellian" # Distribution function (maxwellian, jkp, cgl) +f0type = "maxwellian" # Distribution function (maxwellian, park, cgl) atol_xlmda = 1e-9 # Absolute tolerance for inner pitch + energy integrations rtol_xlmda = 1e-5 # Relative tolerance for inner pitch + energy integrations diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index c92d6e698..8448acbce 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -162,7 +162,7 @@ end Plot the edge stability scan energy components (et, ep, ev, evonly) vs ψ_N. The edge scan evaluates `δW_total = δW_plasma + δW_vacuum` at each stored integration step -in the region [psiedge, psilim], with the plasma boundary swept from psiedge to psilim. +in the region [dW_edge_scan_start, psilim], with the plasma boundary swept from dW_edge_scan_start to psilim. A positive et indicates stability; the truncation point is chosen at the peak et. Four subplots are shown: @@ -177,7 +177,7 @@ A horizontal dashed line at zero marks the stability boundary. A vertical dashed ### Arguments - - `h5path`: Path to a GPEC HDF5 output file produced with `psiedge < psilim` + - `h5path`: Path to a GPEC HDF5 output file produced with `dW_edge_scan_start < psilim` ### Keyword arguments @@ -204,7 +204,7 @@ function plot_edge_stability_scan(h5path; save_path=nothing, ylims=(-2, 3), kwar end if !has_scan - @warn "No edge_scan group in $h5path. Run with psiedge < psilim to generate it." + @warn "No edge_scan group in $h5path. Run with dW_edge_scan_start < psilim to generate it." return nothing end diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index ed6582131..66976648b 100644 --- a/src/Equilibrium/DirectEquilibrium.jl +++ b/src/Equilibrium/DirectEquilibrium.jl @@ -434,7 +434,7 @@ function _build_psi_grid(equil_params, psilow, psihigh) N_core = round(Int, mpsi * log_core / log_total) N_mid = mpsi - N_edge - N_core make_optimal_psi_grid(psilow, psihigh, N_core, N_mid, N_edge) - elseif equil_params.grid_type == "ldp" + elseif equil_params.grid_type == "rational_packed" [psilow + (psihigh - psilow) * sin((ipsi / mpsi) * (π / 2))^2 for ipsi in 0:mpsi] elseif equil_params.grid_type == "pow1" # Fortran powspace(psilow, psihigh, 1, mpsi+1, "upper") — edge-packed grid (equil/grid.f90:92-195) @@ -532,7 +532,7 @@ robustness. rewind!(pool, Float64) end - # Temporary splines for q0 extrapolation and optional newq0 revision + # Temporary splines for q0 extrapolation and optional q0_override revision profiles = ProfileSplines( psi_nodes, sq_nodes[:, 1], # F * 2π @@ -543,17 +543,17 @@ robustness. # q(0) by linear extrapolation from innermost surface q0 = profiles.q_spline.y[1] - profiles.q_deriv(psi_nodes[1]; hint=Ref(1)) * psi_nodes[1] if q0 <= 0.0 - @warn "q0 extrapolation to axis gives q0 = $(@sprintf("%.3f", q0)) ≤ 0 — likely a spline artifact from psilow being too large; check psilow or use newq0 to override." + @warn "q0 extrapolation to axis gives q0 = $(@sprintf("%.3f", q0)) ≤ 0 — likely a spline artifact from psilow being too large; check psilow or use q0_override to override." end - if equil_params.newq0 == -1 - equil_params.newq0 = -q0 + if equil_params.q0_override == -1 + equil_params.q0_override = -q0 end - if equil_params.newq0 != 0.0 - @info "Revising q-profile for newq0 = $(@sprintf("%.3f", equil_params.newq0))" + if equil_params.q0_override != 0.0 + @info "Revising q-profile for q0_override = $(@sprintf("%.3f", equil_params.q0_override))" f0 = profiles.F_spline.y[1] - profiles.F_deriv(psi_nodes[1]; hint=Ref(1)) * psi_nodes[1] - f0fac = f0^2 * ((equil_params.newq0 / q0)^2 - 1.0) + f0fac = f0^2 * ((equil_params.q0_override / q0)^2 - 1.0) for i in 1:(mpsi+1) - ffac = sqrt(1.0 + f0fac / profiles.F_spline.y[i]^2) * sign(equil_params.newq0) + ffac = sqrt(1.0 + f0fac / profiles.F_spline.y[i]^2) * sign(equil_params.q0_override) sq_nodes[i, 1] *= ffac sq_nodes[i, 4] *= ffac rzphi_nodes[i, :, 3] .*= ffac diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl index 4bb29f8f2..e46a30372 100644 --- a/src/Equilibrium/EquilibriumTypes.jl +++ b/src/Equilibrium/EquilibriumTypes.jl @@ -24,8 +24,8 @@ Bundles all necessary settings originally specified in the equil fortran namelis - `r0exp::Float64` - Major radius normalization for CHEASE/EQDSK [m] - `b0exp::Float64` - On-axis toroidal field normalization for CHEASE/EQDSK [T] - `grid_type::String` - Grid type for flux surface discretization ("auto" — two-pass measured-curvature - refinement when mpsi=0, three-region log layout when mpsi>0; "ldp", "pow1", "uniform"; - "log_asymptotic" is a legacy alias for "auto") + refinement when mpsi=0, three-region log layout when mpsi>0; "rational_packed", "pow1", "uniform"; + "log_asymptotic" is a legacy alias for "auto"; "ldp" is a deprecated alias for "rational_packed") - `psilow::Float64` - Lower limit of normalized flux coordinate - `psihigh::Float64` - Upper limit of normalized flux coordinate - `mpsi::Int` - Number of radial grid intervals; 0 with grid_type="auto" selects the @@ -35,10 +35,10 @@ Bundles all necessary settings originally specified in the equil fortran namelis - `psi_accuracy::Float64` - Target relative accuracy τ of splined profile derivatives for the two-pass auto grid (knot count scales as τ^(-1/3)) - `mtheta::Int` - Number of poloidal grid points - - `newq0::Int` - Override for on-axis safety factor (0 = use input value) + - `q0_override::Float64` - Override for the on-axis safety factor q0 (0 = use input value; -1 = flip the sign of the extrapolated q0) - `etol::Float64` - Error tolerance for equilibrium solver - `force_termination::Bool` - Terminate after equilibrium setup (skip stability calculations) - - `use_galgrid::Bool` - Use the same grid as galerkin method + - `use_galerkin_grid::Bool` - Use the same grid as galerkin method """ @kwdef mutable struct EquilibriumConfig eq_type::String = "efit" @@ -64,11 +64,11 @@ Bundles all necessary settings originally specified in the equil fortran namelis psi_accuracy::Float64 = 0.001 mtheta::Int = 512 - newq0::Int = 0 + q0_override::Float64 = 0.0 etol::Float64 = 1e-10 force_termination::Bool = false - use_galgrid::Bool = true + use_galerkin_grid::Bool = true # IMAS-specific: expected COCOS convention of the input dd.equilibrium (11=IMAS standard, 2=GPEC internal) imas_cocos::Int = 11 @@ -81,8 +81,8 @@ Bundles all necessary settings originally specified in the equil fortran namelis # so their incoming values are ignored (hence `_`). function EquilibriumConfig(eq_type, eq_filename, r0exp, b0exp, jac_type, _, _, _, _, jac_custom_power_bp, jac_custom_power_b, jac_custom_power_r, jac_custom_power_rc, - grid_type, psilow, psihigh, mpsi, psi_accuracy, mtheta, newq0, etol, - force_termination, use_galgrid, imas_cocos) + grid_type, psilow, psihigh, mpsi, psi_accuracy, mtheta, q0_override, etol, + force_termination, use_galerkin_grid, imas_cocos) if jac_type == "hamada" @info "Forcing hamada coordinate jacobian exponents: power_*" power_b = 0 @@ -147,8 +147,8 @@ Bundles all necessary settings originally specified in the equil fortran namelis psihigh = min(psihigh, 1.0) return new(eq_type, eq_filename, r0exp, b0exp, jac_type, power_bp, power_b, power_r, power_rc, jac_custom_power_bp, jac_custom_power_b, jac_custom_power_r, jac_custom_power_rc, - grid_type, psilow, psihigh, mpsi, psi_accuracy, mtheta, newq0, etol, - force_termination, use_galgrid, imas_cocos) + grid_type, psilow, psihigh, mpsi, psi_accuracy, mtheta, q0_override, etol, + force_termination, use_galerkin_grid, imas_cocos) end end @@ -176,6 +176,12 @@ function EquilibriumConfig(equil_dict::Dict{String,Any}, base_path::String="./") end end + # Deprecated grid_type value: initials "ldp" renamed to the descriptive spelling. + if get(config_data, "grid_type", "") == "ldp" + @warn "grid_type = \"ldp\" in [Equilibrium] is deprecated; use grid_type = \"rational_packed\". The old spelling will be removed after v2.0.0." + config_data["grid_type"] = "rational_packed" + end + # Construct validated struct config = EquilibriumConfig(; symbolize_keys(config_data)...) # Only resolve `eq_filename` against `base_path` if the user actually diff --git a/src/Equilibrium/InverseEquilibrium.jl b/src/Equilibrium/InverseEquilibrium.jl index 725391b15..a7d69c8d4 100644 --- a/src/Equilibrium/InverseEquilibrium.jl +++ b/src/Equilibrium/InverseEquilibrium.jl @@ -60,7 +60,7 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No mtheta = config.mtheta psilow = config.psilow psihigh = config.psihigh - newq0 = config.newq0 + q0_override = config.q0_override # c----------------------------------------------------------------------- # c allocate and define local arrays. @@ -155,7 +155,7 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No N_core = round(Int, mpsi * log_core / log_total) N_mid = mpsi - N_edge - N_core sq_xs = make_optimal_psi_grid(psilow, psihigh, N_core, N_mid, N_edge) - elseif grid_type == "ldp" + elseif grid_type == "rational_packed" if mpsi == 0 mpsi = 128 end @@ -295,16 +295,16 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No f1_sq_lo = sq_deriv(sq_xs[1]) f1_sq_hi = sq_deriv(sq_xs[end]) q0 = f_sq[1, 4] - f1_sq_lo[4] * sq_xs[1] - if newq0 == -1 - newq0 = -q0 + if q0_override == -1 + q0_override = -q0 end - if newq0 != 0 + if q0_override != 0 f0 = f_sq[1, 2] - f1_sq_lo[2] * sq_xs[1] - f0fac = f0^2 * ((newq0 / q0)^2 - 1) - q0 = newq0 + f0fac = f0^2 * ((q0_override / q0)^2 - 1) + q0 = q0_override for ipsi in 0:mpsi - ffac = sqrt(1 + f0fac / f_sq[ipsi+1, 1]^2) * sign(newq0) + ffac = sqrt(1 + f0fac / f_sq[ipsi+1, 1]^2) * sign(q0_override) sq_fs[ipsi+1, 1] *= ffac sq_fs[ipsi+1, 4] *= ffac rzphi_fs[ipsi+1, :, 3] *= ffac diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index f563dca45..42305b71a 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -84,10 +84,10 @@ function balance_integration_chunks(chunks::Vector{IntegrationChunk}, ctrl::Forc # assemble_fm_matrix(condition=true) can't keep accumulated products well-conditioned # because single long-span propagators may already have cond ~ 10²⁴. min_bvp_intervals = 8 * (intr.msing + 1) + intr.msing - # Use the effective parallel width (capped by ctrl.parallel_threads) rather than - # Threads.nthreads() — otherwise a user on `julia -t 16` who sets parallel_threads=2 + # Use the effective parallel width (capped by ctrl.integrator_threads) rather than + # Threads.nthreads() — otherwise a user on `julia -t 16` who sets integrator_threads=2 # for determinism still pays for 4× the requested sub-chunk count. - effective_threads = min(Threads.nthreads(), max(ctrl.parallel_threads, 1)) + effective_threads = min(Threads.nthreads(), max(ctrl.integrator_threads, 1)) target_n = max(min_chunks, 4 * effective_threads, min_bvp_intervals) result = collect(chunks) @@ -140,21 +140,21 @@ end eulerlagrange_integration(ctrl, equil, ffit, intr) -> (odet, propagators, chunks, S_left) Integrate the Euler-Lagrange equations from the axis to `intr.psilim`, crossing each singular -surface on the way (Fortran `ode_run`). Dispatches on `ctrl` to the parallel propagator BVP -(`use_parallel`), the dual Riccati formulation (`use_riccati`), or -[`serial_eulerlagrange_integration`](@ref). +surface on the way (Fortran `ode_run`). Dispatches on `ctrl.integrator`: the STRIDE propagator +BVP (`"stride"`), the dual Riccati formulation (`"riccati"`), or +[`serial_eulerlagrange_integration`](@ref) (`"serial"`). -Only the parallel branch populates `propagators` / `chunks` / `S_left`, which +Only the stride branch populates `propagators` / `chunks` / `S_left`, which `compute_delta_prime_matrix!` consumes for the Δ' BVP; the other two return `nothing` for all three. """ function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) - # Dispatch to parallel or Riccati solver if requested. - # Parallel path returns (odet, propagators, chunks, S_at_surface_left) for deferred Δ' BVP. - if ctrl.use_parallel + # Dispatch on the integrator algorithm. + # Stride path returns (odet, propagators, chunks, S_at_surface_left) for deferred Δ' BVP. + if ctrl.integrator == "stride" return parallel_eulerlagrange_integration(ctrl, equil, ffit, intr) - elseif ctrl.use_riccati + elseif ctrl.integrator == "riccati" return (riccati_eulerlagrange_integration(ctrl, equil, ffit, intr), nothing, nothing, nothing) end return serial_eulerlagrange_integration(ctrl, equil, ffit, intr) @@ -166,7 +166,7 @@ end Serial shooting branch of [`eulerlagrange_integration`](@ref): integrates chunk by chunk, applying Gaussian reduction whenever a solution norm ratio exceeds `ctrl.ucrit` and undoing it via `transform_u!` at the end, so `odet.u_store` comes back dense in the axis basis. Call -directly to force this branch regardless of `ctrl.use_parallel` / `ctrl.use_riccati`; `verbose` +directly to force this branch regardless of `ctrl.integrator`; `verbose` overrides `ctrl.verbose` for progress logging. """ function serial_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal; @@ -214,7 +214,7 @@ function serial_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::E odet.step -= 1 trim_storage!(odet) - # Edge-dW scan over [psiedge, psilim] — populates odet.edge_scan for HDF5 output. + # Edge-dW scan over [dW_edge_scan_start, psilim] — populates odet.edge_scan for HDF5 output. # The scan mutates odet.psifac and odet.u internally; save/restore them around the call. # findmax_dW_edge! also (re)allocates odet.edge_scan; that field is the diagnostic # product and is intentionally NOT restored. @@ -224,7 +224,7 @@ function serial_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::E # location. Legacy path (true) reproduces the ode_record_edge heuristic from Fortran # STRIDE — psilim/qlim/u are pulled back to the dW peak. Preserved for experimental # work; see docstring in ForceFreeStatesStructs.jl for the reliability caveats. - if ctrl.psiedge < intr.psilim + if ctrl.dW_edge_scan_start < intr.psilim saved_psifac, saved_u = odet.psifac, copy(odet.u) peak_step = findmax_dW_edge!(odet, ctrl, equil, ffit, intr) if ctrl.truncate_at_dW_peak @@ -752,7 +752,7 @@ function integrate_el_region!( near_start = abs(odet.q - q_start) < near_q_frac * q_range || steps_in_segment[] == 1 near_end = abs(odet.q - q_end) < near_q_frac * q_range # Always save in the edge scan region so findmax_dW_edge! has dense q coverage. - in_edge_scan = ctrl.psiedge < intr.psilim && integrator.t >= ctrl.psiedge + in_edge_scan = ctrl.dW_edge_scan_start < intr.psilim && integrator.t >= ctrl.dW_edge_scan_start if near_start || near_end || (odet.total_steps % ctrl.save_interval == 0) || in_edge_scan # q at the accepted point, not the last internal Runge-Kutta stage @@ -881,7 +881,7 @@ end """ findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) -Records the total dW in the integration region between `ctrl.psiedge` and +Records the total dW in the integration region between `ctrl.dW_edge_scan_start` and `ctrl.psilim`. This performs the same function as `ode_record_edge` in the Fortran, but everything is now done post-integration which cleans up the logic, i.e. no "_edge" arrays. @@ -897,8 +897,8 @@ for clarity. We create the wv matrix spline once prior to the loop. """ function findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) - # Find the first ODE step at or past psiedge; all subsequent steps are contiguous edge steps - edge_start = findfirst(i -> odet.psi_store[i] >= ctrl.psiedge, 1:odet.step) + # Find the first ODE step at or past dW_edge_scan_start; all subsequent steps are contiguous edge steps + edge_start = findfirst(i -> odet.psi_store[i] >= ctrl.dW_edge_scan_start, 1:odet.step) N_edge = odet.step - edge_start + 1 # Initialize EdgeScanState sized exactly to the number of edge steps @@ -908,7 +908,7 @@ function findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::E es.psi .= odet.psi_store[edge_start:odet.step] es.q .= odet.q_store[edge_start:odet.step] - # Create a rough spline for wv matrix between psiedge -> psilim so we can approximate dW + # Create a rough spline for wv matrix between dW_edge_scan_start -> psilim so we can approximate dW es.wvmat = free_compute_wv_spline(ctrl, equil, intr) # Loop with compact index j into EdgeScanState; ODE index is edge_start + j - 1. diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index d9c4770d2..93e79a56f 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -206,7 +206,7 @@ A mutable struct holding internal state variables for stability calculations. via `pest3_decompose(dp_raw)` — needed for the full det(D' − D(γ)) = 0 eigenvalue problem with Glasser stabilization. - Empty unless `ctrl.use_parallel` is true. No ½ prefactor is applied (matches + Empty unless `ctrl.integrator == "stride"`. No ½ prefactor is applied (matches Fortran rdcon; Pletzer–Dewar paper multiplies by ½). """ delta_prime_raw::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) @@ -230,7 +230,6 @@ gpec.toml. - `nn_high::Int` - Upper bound for toroidal modes - `delta_mlow::Int` - Expands lower bound of Fourier harmonics by delta_mlow - `delta_mhigh::Int` - Expands upper bound of Fourier harmonics by delta_mhigh - - `nstep::Int` - Maximum number of integration steps (not yet implemented) - `ksing::Int` - Singular surface handling parameter - `eulerlagrange_tolerance::Float64` - Relative tolerance for ODE integration of Euler-Lagrange equations - `ucrit::Float64` - Critical value of unorm ratio to trigger solution normalization. In the standard path it triggers Gaussian reduction; in the Riccati path it triggers `renormalize_riccati_inplace!`. Default `1e4` empirically keeps max(|U₁|, |U₂|) in O(1)–O(10⁴) over the integration domain on DIII-D / Solovev sweeps; lower triggers excess renorms without accuracy gain, higher risks overflow before the next renorm. @@ -245,18 +244,16 @@ gpec.toml. - `kinetic_factor::Float64` - Dimensionless scaling factor for kinetic matrices. Zero (the default) disables the kinetic path; any positive value enables it and scales the kinetic matrices: when kinetic_source="fixed", scales X-shaped test matrices relative to ideal matrix norms; when kinetic_source="calculated", applied as uniform post-hoc multiplier to W and T components. - `qlow::Float64` - Integration terminated at q limit determined by minimum of qlow and q0 from equil - `reform_eq_with_psilim::Bool` - Reform equilibrium with computed psilim (not yet implemented) - - `psiedge::Float64` - If less than psilim, records a dW(ψ) diagnostic scan over [psiedge, psilim] on odet.edge_scan. The integration domain (psilim) is always controlled by qhigh / psihigh and is not modified by this scan (unless `truncate_at_dW_peak=true`, see caveats below). - - `truncate_at_dW_peak::Bool` - When `true` and `psiedge < psilim`, the edge-dW scan's peak location is adopted as the new physical plasma edge — `intr.psilim`/`intr.qlim`/`odet.u` are pulled back to the peak, AND the FM Δ' chunks/propagators are made self-consistent with the new boundary (the chunk that straddles the peak is rebuilt + re-integrated; any chunks past the peak are dropped). This reproduces the spirit of the original ode_record_edge heuristic from Fortran STRIDE while keeping Δ' and δW well-defined at the new boundary. The Δ' metric is still physically dependent on where the peak falls in the edge band, so use this flag deliberately when you mean to scan against the peak-defined edge (e.g. for studying edge-mode regimes); leave at `false` (default) for the full-domain Δ' at `qhigh` / `psihigh` / `dmlim`. + - `dW_edge_scan_start::Float64` - If less than psilim, records a dW(ψ) diagnostic scan over [dW_edge_scan_start, psilim] on odet.edge_scan. The integration domain (psilim) is always controlled by qhigh / psihigh and is not modified by this scan (unless `truncate_at_dW_peak=true`, see caveats below). + - `truncate_at_dW_peak::Bool` - When `true` and `dW_edge_scan_start < psilim`, the edge-dW scan's peak location is adopted as the new physical plasma edge — `intr.psilim`/`intr.qlim`/`odet.u` are pulled back to the peak, AND the FM Δ' chunks/propagators are made self-consistent with the new boundary (the chunk that straddles the peak is rebuilt + re-integrated; any chunks past the peak are dropped). This reproduces the spirit of the original ode_record_edge heuristic from Fortran STRIDE while keeping Δ' and δW well-defined at the new boundary. The Δ' metric is still physically dependent on where the peak falls in the edge band, so use this flag deliberately when you mean to scan against the peak-defined edge (e.g. for studying edge-mode regimes); leave at `false` (default) for the full-domain Δ' at `qhigh` / `psihigh` / `dmlim`. - `diagnose::Bool` - Enable diagnostic output (not yet implemented) - - `diagnose_ca::Bool` - Enable asymptotic coefficient diagnostics (not yet implemented) - `write_outputs_to_HDF5::Bool` - Write results to HDF5 format - `HDF5_filename::String` - Name of HDF5 output file - `save_interval::Int` - Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. (Same as `euler_step` in the Fortran) - `force_termination::Bool` - Terminate after force-free states (skip perturbed equilibrium calculations) - - `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 at the end of the propagator BVP and let it replace the `odet` returned to the main pipeline. This populates `u_store` / `du_store` / `xi_s_store` densely in the axis (EL) basis — the only convention the PerturbedEquilibrium / FieldReconstruction downstream code consumes correctly. Without it the parallel path stores only chunk-endpoint Riccati S matrices with diagnostic derivatives (see Riccati.jl docstring caveats), and HDF5 `integration/xi_psi`/`dxi_psi`/`xi_s` are unusable. Δ' (`singular/delta_prime_matrix`) is computed from the parallel BVP and is bit-identical between `populate_dense_xi=true` and `false`. Energies (`vacuum/ep`/`ev`/`et`) are computed by `free_run` from `odet`, so with `populate_dense_xi=true` they match what a pure serial run (`use_parallel=false`) would produce; with `populate_dense_xi=false` they use the parallel-pass Riccati `odet.u` instead (differs by the ~0.12 % Riccati-vs-axis algorithmic gap on DIIID-class cases). **Default `false`** to avoid paying the dense-pass cost on Δ'/vacuum/ideal-stability-only runs; **PerturbedEquilibrium-using configs must set `populate_dense_xi = true` explicitly** when `use_parallel = true` (otherwise PE silently reads Riccati-basis garbage). Auto-disabled when `force_termination = true` regardless of the user setting, since the dense pass has no downstream consumer in that case. Approximate cost when enabled: one extra serial EL integration (~1× the parallel BVP wall-clock for typical N). + - `integrator::String` - Which Euler-Lagrange integration algorithm to use. `"stride"` (default): fundamental-matrix (propagator) chunk integration following Fortran STRIDE, assembled via a BVP — the only path that produces the Δ' matrix (`SingularSurfaces/delta_prime_matrix`) consumed by SLAYER/GGJ downstream. Each chunk is integrated independently from identity IC using `Threads.@threads` (thread cap `integrator_threads`; `integrator_threads = 1` runs the chunks serially and bit-deterministically). Requires `singfac_min != 0`; crossings use the Riccati-style algorithm (no Gaussian reduction). `"riccati"`: 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. `"serial"`: the serial shooting method with Gaussian reduction (see `serial_eulerlagrange_integration`). + - `integrator_threads::Int` - Cap on the number of threads the stride-integrator 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 `integrator_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 stride path produces bit-identical Δ′ across thread counts; `integrator_threads = 2` is about 20% faster than serial and saturates the speedup. If a stride run diverges, drop to `integrator_threads = 1` rather than switching to `integrator = "serial"` — the latter is silently wrong. Capped at `Threads.nthreads()`. + - `populate_dense_xi::Bool` - When `integrator = "stride"`, append a serial Euler-Lagrange pass at the end of the propagator BVP and let it replace the `odet` returned to the main pipeline. This populates `u_store` / `du_store` / `xi_s_store` densely in the axis (EL) basis — the only convention the PerturbedEquilibrium / FieldReconstruction downstream code consumes correctly. Without it the stride path stores only chunk-endpoint Riccati S matrices with diagnostic derivatives (see Riccati.jl docstring caveats), and the HDF5 `ForceFreeStates/Solutions/ForwardIntegration` `xi_psi`/`dxi_psi`/`xi_s` are unusable. Δ' (`SingularSurfaces/delta_prime_matrix`) is computed from the stride BVP and is bit-identical between `populate_dense_xi=true` and `false`. Energies (`ForceFreeStates/FreeBoundaryStability` `eigenmode_*_energies`) are computed by `free_run` from `odet`, so with `populate_dense_xi=true` they match what a pure serial run (`integrator = "serial"`) would produce; with `populate_dense_xi=false` they use the stride-pass Riccati `odet.u` instead (differs by the ~0.12 % Riccati-vs-axis algorithmic gap on DIIID-class cases). **Default `false`** to avoid paying the dense-pass cost on Δ'/vacuum/ideal-stability-only runs; **PerturbedEquilibrium-using configs must set `populate_dense_xi = true` explicitly** when `integrator = "stride"` (otherwise PE silently reads Riccati-basis garbage). Auto-disabled when `force_termination = true` regardless of the user setting, since the dense pass has no downstream consumer in that case. Approximate cost when enabled: one extra serial EL integration (~1× the stride BVP wall-clock for typical N). - `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 @@ -270,13 +267,12 @@ gpec.toml. nn_high::Int = 0 delta_mlow::Int = 0 delta_mhigh::Int = 0 - nstep::Int = typemax(Int) ksing::Int = -1 eulerlagrange_tolerance::Float64 = 1e-8 ucrit::Float64 = 1e4 numsteps_init::Int = 4000 numunorms_init::Int = 100 - singfac_min::Float64 = 1e-4 # Matches Fortran STRIDE; required nonzero for use_parallel path. + singfac_min::Float64 = 1e-4 # Matches Fortran STRIDE; required nonzero for the stride integrator path. set_psilim_via_dmlim::Bool = true # Safe default for diverted equilibria (most production use); set false for limited/analytical (LAR, Solovev). Auto-skipped for multi-n. See docstring. dmlim::Float64 = 0.2 sing_order::Int = 6 @@ -285,18 +281,16 @@ gpec.toml. kinetic_factor::Float64 = 0.0 qlow::Float64 = 0.0 reform_eq_with_psilim::Bool = false - psiedge::Float64 = 0.99 + dW_edge_scan_start::Float64 = 0.99 truncate_at_dW_peak::Bool = false # Edge-dW peak becomes new physical edge; Δ' BVP made self-consistent. See docstring. - parallel_threads::Int = 2 + integrator_threads::Int = 2 diagnose::Bool = false - diagnose_ca::Bool = false write_outputs_to_HDF5::Bool = true HDF5_filename::String = "gpec.h5" save_interval::Int = 3 force_termination::Bool = false - use_riccati::Bool = false - use_parallel::Bool = true # Default on: unlocks singular/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). + integrator::String = "stride" # Default stride: unlocks SingularSurfaces/delta_prime_matrix (STRIDE BVP Δ' matrix) used by SLAYER/GGJ downstream. + populate_dense_xi::Bool = false # When integrator="stride", set to true ONLY if a PerturbedEquilibrium pipeline will consume dense ξ. Default false avoids the ~1× stride-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). # --- RDCON outer-region Galerkin Δ′ solver (gal_solve port) --- @@ -447,7 +441,7 @@ end """ EdgeScanState -Holds the state and results for the edge dW stability scan over ψ ∈ [psiedge, psilim]. +Holds the state and results for the edge dW stability scan over ψ ∈ [dW_edge_scan_start, psilim]. Initialized and populated by `findmax_dW_edge!`; results written to HDF5 under `EdgeScan/`. The energies are generalized (W, N) pencil values: power-normalized and invariant to the working (Jacobian) coordinate (see `power_norm_matrix!`). @@ -562,7 +556,7 @@ and a small set of temporary matrices and factors used to compute singular-layer - `index::Array{Int,2}` - Index matrix used for sorting solution norms with shape `(numpert_total, numunorms_init)`. - - `sing_flag::Vector{Bool}` - Boolean flags indicating which stored normalizations correspond to singular solutions # Edge dW scan state and results (disabled sentinel when psiedge >= psilim, i.e. no edge scan) + - `sing_flag::Vector{Bool}` - Boolean flags indicating which stored normalizations correspond to singular solutions # Edge dW scan state and results (disabled sentinel when dW_edge_scan_start >= psilim, i.e. no edge scan) (length `numunorms_init`). - `zeroed_idx::Vector{Vector{Int}}` - For each ideal rational surface jump, a vector of indices of solutions that were zeroed. # Data for integrator @@ -593,7 +587,7 @@ and a small set of temporary matrices and factors used to compute singular-layer ca_r::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, msing) ca_l::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, msing) - # Edge dW scan state and results (disabled sentinel when psiedge >= psilim, i.e. no edge scan) + # Edge dW scan state and results (disabled sentinel when dW_edge_scan_start >= psilim, i.e. no edge scan) edge_scan::EdgeScanState = EdgeScanState(numpert_total, 0) # Data for integrator diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index adf545b1c..e7ed75f68 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -170,7 +170,7 @@ q-window minimum. # Number of psi grid points for the spline: 4 per q-window minimum # TODO: 4 spline points is arbitrary - is there a better way? - qedge = profiles.q_spline(ctrl.psiedge) + qedge = profiles.q_spline(ctrl.dW_edge_scan_start) npsi = max(4, ceil(Int, (intr.qlim - qedge) * intr.nhigh * 4)) psi_array = zeros!(pool, Float64, npsi + 1) wv_array = zeros!(pool, ComplexF64, npsi + 1, intr.numpert_total, intr.numpert_total) @@ -179,7 +179,7 @@ q-window minimum. # Space points evenly in q over [qedge, qlim] (i=1 → qedge, i=npsi+1 → qlim) qi = qedge + (intr.qlim - qedge) * ((i - 1) / npsi) - psii = ctrl.psiedge + (intr.psilim - ctrl.psiedge) * ((i - 1) / npsi) + psii = ctrl.dW_edge_scan_start + (intr.psilim - ctrl.dW_edge_scan_start) * ((i - 1) / npsi) psi_array[i] = find_zero( (psi -> profiles.q_spline(psi) - qi, psi -> profiles.q_deriv(psi)), @@ -207,7 +207,7 @@ end Compute total complex energy eigenvalue (total1). This is a trimmed down version of `free_run` that only computes the total energy eigenvalue for the mode unstable mode, used in `findmax_dW_edge!` -which calls this function at each step in the psiedge -> psilim region of integration. This performs +which calls this function at each step in the dW_edge_scan_start -> psilim region of integration. This performs the same function as `free_test` in the Fortran code, except we have moved the creation of the wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvmat` (a complex-valued spline). """ diff --git a/src/ForceFreeStates/Riccati.jl b/src/ForceFreeStates/Riccati.jl index 9e86aadcf..68f658db1 100644 --- a/src/ForceFreeStates/Riccati.jl +++ b/src/ForceFreeStates/Riccati.jl @@ -1307,7 +1307,7 @@ otherwise mirrors, differing in three places: renormalizes to (S_new, I) in one step 3. `transform_u!` is skipped — S is already the true solution, so there is no reduction to undo -Enable via `use_riccati = true` in the `[ForceFreeStates]` section of gpec.toml. +Enable via `integrator = "riccati"` in the `[ForceFreeStates]` section of gpec.toml. """ function riccati_eulerlagrange_integration( ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, @@ -1353,13 +1353,13 @@ function riccati_eulerlagrange_integration( end end - # Edge-dW scan over [psiedge, psilim] — populates odet.edge_scan for HDF5 output. + # Edge-dW scan over [dW_edge_scan_start, psilim] — populates odet.edge_scan for HDF5 output. # See EulerLagrange.jl counterpart and ForceFreeStatesControl docstring for the # diagnostic vs legacy-truncation semantics and reliability caveats on # truncate_at_dW_peak=true. odet.step -= 1 trim_storage!(odet) - if ctrl.psiedge < intr.psilim + if ctrl.dW_edge_scan_start < intr.psilim saved_psifac, saved_u = odet.psifac, copy(odet.u) peak_step = findmax_dW_edge!(odet, ctrl, equil, ffit, intr) if ctrl.truncate_at_dW_peak @@ -1622,7 +1622,7 @@ concurrently using `Threads.@threads`, then re-integrates the outer plasma seria without renormalization); Riccati integration keeps matrices bounded and provides dense checkpoints for `findmax_dW_edge!`. -Enable via `use_parallel = true` in `[ForceFreeStates]` of gpec.toml. Requires `singfac_min != 0`. +Enable via `integrator = "stride"` in `[ForceFreeStates]` of gpec.toml. Requires `singfac_min != 0`. **Key differences from serial integration:** - No Gaussian reduction in the propagator BVP phase (crossings use the @@ -1650,7 +1650,7 @@ function parallel_eulerlagrange_integration( ) odet = _initialize_parallel_odet(ctrl, equil, ffit, intr) chunks, propagators, odet_proxies = _setup_parallel_chunks_and_proxies(odet, ctrl, intr) - bvp_threads = max(1, min(Threads.nthreads(), ctrl.parallel_threads)) + bvp_threads = max(1, min(Threads.nthreads(), ctrl.integrator_threads)) _log_parallel_start(ctrl, odet, equil, chunks, bvp_threads) _run_parallel_bvp_phase!(propagators, chunks, ctrl, equil, ffit, intr, odet_proxies, bvp_threads) @@ -1725,14 +1725,14 @@ function _log_parallel_start(ctrl::ForceFreeStatesControl, odet::OdeState, chunks::Vector{IntegrationChunk}, bvp_threads::Int) ctrl.verbose || return @info " ψ = $((@sprintf "%.3f" odet.psifac)), q = $((@sprintf "%.3f" equil.profiles.q_spline(odet.psifac)))" - @info " Parallel FM: $(length(chunks)) chunks, $bvp_threads BVP thread$(bvp_threads == 1 ? "" : "s") (julia_nthreads=$(Threads.nthreads()), ctrl.parallel_threads=$(ctrl.parallel_threads))" + @info " Parallel FM: $(length(chunks)) chunks, $bvp_threads BVP thread$(bvp_threads == 1 ? "" : "s") (julia_nthreads=$(Threads.nthreads()), ctrl.integrator_threads=$(ctrl.integrator_threads))" end # Integrate each chunk's FM propagator from identity IC. Serial when bvp_threads == 1 # (bit-deterministic; ~20% slower than 2-thread but immune to thread- # schedule sensitivity). Parallel uses :static scheduler so Threads.threadid() returns a # stable index into odet_proxies. If a parallel run ever diverges on a delicate equilibrium, -# drop to parallel_threads = 1 rather than use_parallel = false — the latter is silently wrong. +# drop to integrator_threads = 1 rather than integrator = "serial" — the latter is silently wrong. function _run_parallel_bvp_phase!(propagators::Vector{ChunkPropagator}, chunks::Vector{IntegrationChunk}, ctrl::ForceFreeStatesControl, @@ -1827,7 +1827,7 @@ function _reintegrate_outer_plasma!(odet::OdeState, last_crossing_step::Int, # Post: odet.u is in (S, I) form; odet.step points to next empty slot. end -# Edge-dW scan over [psiedge, psilim] — populates odet.edge_scan for HDF5. By default +# Edge-dW scan over [dW_edge_scan_start, psilim] — populates odet.edge_scan for HDF5. By default # (truncate_at_dW_peak=false) it's diagnostic-only: integration domain is unchanged. # When truncate_at_dW_peak=true, the dW peak becomes the new physical edge: intr.psilim, # odet, propagators, and chunks are made self-consistent (straddling chunk rebuilt with @@ -1843,7 +1843,7 @@ function _handle_edge_dW_scan!(odet::OdeState, chunks::Vector{IntegrationChunk}, N = intr.numpert_total odet.step -= 1 trim_storage!(odet) - ctrl.psiedge < intr.psilim || return chunks, propagators + ctrl.dW_edge_scan_start < intr.psilim || return chunks, propagators saved_psifac, saved_u = odet.psifac, copy(odet.u) peak_step = findmax_dW_edge!(odet, ctrl, equil, ffit, intr) @@ -1924,7 +1924,7 @@ which `compute_delta_prime_matrix!` uses). Called from `parallel_eulerlagrange_integration` when `ctrl.populate_dense_xi = true`. Approximate cost: one serial EL integration on top of the parallel BVP phase. Required to make -`use_parallel = true` produce DCON eigenfunctions usable by the +`integrator = "stride"` produce DCON eigenfunctions usable by the PerturbedEquilibrium downstream pipeline. """ function _populate_dense_xi_via_serial_el!( diff --git a/src/ForcingTerms/CoilFourier.jl b/src/ForcingTerms/CoilFourier.jl index 62996bf59..dbc1a160c 100644 --- a/src/ForcingTerms/CoilFourier.jl +++ b/src/ForcingTerms/CoilFourier.jl @@ -62,7 +62,7 @@ Defaults to the **outermost computed surface** `equil.rzphi_xs[end]` (i.e. `psih square root: `rfac = SQRT(crzphi_f(1))` (`coil/field.F:170`, where `crzphi_f(1)` is `rzphi_rsquared` and `crzphi_f(2)` is `rzphi_offset`); `coil/field.F:133` calls that mesh the "control surface mesh". `psilim = psihigh` (`dcon/sing.f:170`) and is only ever moved *inward* -by `sas_flag`/`qhigh`/`psiedge` truncation. Since `psihigh` is the last knot of the radial +by `sas_flag`/`qhigh`/`dW_edge_scan_start` truncation. Since `psihigh` is the last knot of the radial grid these splines are built on (`equil/inverse.f:142`), Fortran evaluates exactly ON the last knot and never extrapolates. The docs state it directly: the external field is specified "on the surface of the GPEC plasma boundary defined by the psihigh variable in equil.in" @@ -83,7 +83,7 @@ auto grid was refined. NOTE ON THE DEFAULT: the physically correct control surface is `psilim`, the *integration* limit, not `psihigh`, the *equilibrium spline* limit. They are equal unless -`dmlim`/`qhigh`/`psiedge` truncation fires, in which case `psilim < psihigh`. PPPL shipped a fix +`dmlim`/`qhigh`/`dW_edge_scan_start` truncation fires, in which case `psilim < psihigh`. PPPL shipped a fix for exactly this confusion (`docs/releases.rst:281`: "Fixes inappropriate uses of psihigh, which may not be the end of integration psilim if sas_flag, qhigh, or peak_flag are used"). diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 197a4c33f..18306b73c 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -80,9 +80,14 @@ using .ForceFreeStates: find_kinetic_singular_surfaces! using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunctions! using .ForceFreeStates: galerkin_solve, write_galerkin!, GalerkinResult, gal_matched_odestate -const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag") +const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag", "nstep", "diagnose_ca") const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc") +# Old→new TOML key spellings, accepted with a warning until removal after v2.0.0. +const _RENAMED_EQUIL_KEYS = ("newq0" => "q0_override", "use_galgrid" => "use_galerkin_grid") +const _RENAMED_FFS_KEYS = ("parallel_threads" => "integrator_threads", "psiedge" => "dW_edge_scan_start") +const _RENAMED_PE_KEYS = ("reg_spot" => "regularization_width",) + # Drop deprecated keys from a parsed gpec.toml section so legacy files keep parsing # instead of throwing an unknown-keyword error; warn so the removal is not silent. function _drop_deprecated_keys!(table, deprecated_keys, section::String) @@ -95,6 +100,44 @@ function _drop_deprecated_keys!(table, deprecated_keys, section::String) return table end +# Remap old→new key spellings in a parsed gpec.toml section so legacy decks keep working; +# warn on each hit, and let an explicitly set new key win over its old alias. +function _rename_keys!(table, renames, section::String) + for (old, new) in renames + haskey(table, old) || continue + @warn "`$old` in [$section] was renamed to `$new`; the old key is deprecated and will be removed after v2.0.0." + haskey(table, new) || (table[new] = table[old]) + delete!(table, old) + end + return table +end + +# Remap a deprecated value of an enum-like key (old → new spelling), warning on each hit. +function _rename_value!(table, key::String, old, new, section::String) + get(table, key, nothing) == old || return table + @warn "`$key = \"$old\"` in [$section] is deprecated; use `$key = \"$new\"`. The old value will be removed after v2.0.0." + table[key] = new + return table +end + +# The use_parallel/use_riccati boolean pair was replaced by the `integrator` enum; map the +# old flags onto the equivalent algorithm, mirroring the old dispatch order (parallel wins +# over riccati, both false means serial) including the old use_parallel=true default. +function _remap_integrator_keys!(table) + (haskey(table, "use_parallel") || haskey(table, "use_riccati")) || return table + implied = get(table, "use_parallel", true) ? "stride" : + (get(table, "use_riccati", false) ? "riccati" : "serial") + delete!(table, "use_parallel") + delete!(table, "use_riccati") + if haskey(table, "integrator") + @warn "`use_parallel`/`use_riccati` in [ForceFreeStates] are deprecated and ignored because `integrator` is also set." + else + @warn "`use_parallel`/`use_riccati` in [ForceFreeStates] were replaced by the `integrator` enum; assuming `integrator = \"$implied\"`. The old keys will be removed after v2.0.0." + table["integrator"] = implied + end + return table +end + function main(args::Vector{String}=String[]; dd::Union{IMASdd.dd,Nothing}=nothing) # Every input source builds a ready `(inputs, eq_config, additional_input)` and hands it to # `main_from_inputs`: a gpec.toml working directory, an IMAS `dd`, or a gpec.h5 snapshot. @@ -133,6 +176,7 @@ function build_inputs_from_toml(path::String; dd::Union{IMASdd.dd,Nothing}=nothi inputs = TOML.parsefile(joinpath(path, "gpec.toml")) haskey(inputs, "Equilibrium") || error("No [Equilibrium] section in gpec.toml") + _rename_keys!(inputs["Equilibrium"], _RENAMED_EQUIL_KEYS, "Equilibrium") _drop_deprecated_keys!(inputs["Equilibrium"], _DEPRECATED_EQUIL_KEYS, "Equilibrium") eq_config = Equilibrium.EquilibriumConfig(inputs["Equilibrium"], path) @@ -191,8 +235,16 @@ function main_from_inputs( # Build data structures from inputs intr = ForceFreeStatesInternal(; dir_path=path) ffs_table = inputs["ForceFreeStates"] + _rename_keys!(ffs_table, _RENAMED_FFS_KEYS, "ForceFreeStates") + _remap_integrator_keys!(ffs_table) _drop_deprecated_keys!(ffs_table, _DEPRECATED_FFS_KEYS, "ForceFreeStates") ctrl = ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in ffs_table)...) + ctrl.integrator in ("stride", "riccati", "serial") || + error("[ForceFreeStates] integrator = \"$(ctrl.integrator)\" is not one of \"stride\", \"riccati\", \"serial\"") + # SLAYER consumes the STRIDE Δ' matrix; fail at config time rather than after a long run. + if haskey(inputs, "SLAYER") && get(inputs["SLAYER"], "enabled", false) === true && ctrl.integrator != "stride" + error("[SLAYER] requires the STRIDE Δ' matrix: set integrator = \"stride\" in [ForceFreeStates] (got \"$(ctrl.integrator)\")") + end # Determine toroidal mode numbers (n >= 1 required; 0 means "not specified") intr.nlow, intr.nhigh = ctrl.nn_low, ctrl.nn_high @@ -225,6 +277,8 @@ function main_from_inputs( # does not need kinetic_profiles, but the post-PE block always does, so we load # whenever a [KineticForces] section is present or the stability path requests the # calculated source. psio is invariant across grid re-formation. + haskey(inputs, "KineticForces") && + _rename_value!(inputs["KineticForces"], "f0type", "jkp", "park", "KineticForces") kf_ctrl = haskey(inputs, "KineticForces") ? KineticForces.KineticForcesControl(; @@ -582,6 +636,7 @@ function main_from_inputs( ft_ctrl = ForcingTerms.ForcingTermsControl() # Use defaults end + _rename_keys!(inputs["PerturbedEquilibrium"], _RENAMED_PE_KEYS, "PerturbedEquilibrium") pe_ctrl = PerturbedEquilibrium.PerturbedEquilibriumControl(; (Symbol(k) => v for (k, v) in inputs["PerturbedEquilibrium"])... ) @@ -828,7 +883,7 @@ function write_outputs_to_HDF5( out_h5["$fwd/xi_s"] = odet.xi_s_store out_h5["$fwd/crit"] = odet.crit_store - # Write edge stability scan data (only present when psiedge < psilim). + # Write edge stability scan data (only present when dW_edge_scan_start < psilim). # Generalized (W, N) pencil energies — power-normalized, Jacobian-invariant; these are # the values findmax_dW_edge! uses to choose the truncation point. if !isempty(odet.edge_scan.psi) diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 96177903d..ee988238f 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -92,11 +92,11 @@ profiles, without an intermediate file round-trip. callable of `psi` (default `1.0`). - `chi_tor` -- toroidal heat diffusivity [m²/s]. Scalar or a callable of `psi` (default `1.0`). - - `dr_val` -- resistive interchange index `D_R = E + F + H²` + - `delta_crit_D_R` -- resistive interchange index `D_R = E + F + H²` (Glasser-Greene-Johnson 1975) feeding the critical-Δ formulas - (`:lar`, `:rfitzp`, `:toroidal`). When `nothing` (default), Julia + (`:lar`, `:fitzpatrick`, `:toroidal`). When `nothing` (default), Julia derives it per-surface from the equilibrium as - `dr_val_k = D_R(ψ_k) = E_k + F_k + H_k²`, + `delta_crit_D_R_k = D_R(ψ_k) = E_k + F_k + H_k²`, consistent with Connor-Hastie-Helander 2015 (PPCF 57 065001) Eq. 59 which uses `(−D_R)` in the χ_‖-matching critical-Δ. Pass a scalar / vector / callable to override. @@ -106,14 +106,14 @@ profiles, without an intermediate file round-trip. NOT the Mercier index `D_I = E + F + H − 1/4`. The two differ by `(H − 1/2)²`, which is non-trivial on shaped equilibria (~factor 3 on DIII-D); this code uses the physically correct `D_R`. - - `dgeo_val` -- Connor 2015 (PPCF 57 065001) Eq. 59 geometric factor - used by `dc_type=:toroidal`. When `nothing` (default), an error is - raised if `dc_type=:toroidal` is also requested — the auto-derived + - `delta_crit_geo_factor` -- Connor 2015 (PPCF 57 065001) Eq. 59 geometric factor + used by `delta_crit_type=:toroidal`. When `nothing` (default), an error is + raised if `delta_crit_type=:toroidal` is also requested — the auto-derived formula additionally needs ⟨|∇ψ|²⟩ FSA which `ResistGeometry` doesn't currently expose. Pass a scalar / vector / callable to use - a prescribed value. (For `dc_type=:rfitzp` and `:lar`, dgeo_val is + a prescribed value. (For `delta_crit_type=:fitzpatrick` and `:lar`, delta_crit_geo_factor is not consulted.) - - `dc_type` -- `:none` (default), `:lar`, `:rfitzp`, or `:toroidal`. + - `delta_crit_type` -- `:none` (default), `:lar`, `:fitzpatrick`, or `:toroidal`. - `theta` -- poloidal angle at which to measure minor radius (default `0.0`, outboard midplane). - `resistivity_model` -- `SauterNeoModel()` (default), `RedlNeoModel()`, @@ -135,9 +135,9 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; z_i::Real=1.0, chi_perp=1.0, chi_tor=1.0, - dr_val=nothing, - dgeo_val=nothing, - dc_type::Symbol=:none, + delta_crit_D_R=nothing, + delta_crit_geo_factor=nothing, + delta_crit_type::Symbol=:none, theta::Real=0.0, compute_omega_star::Bool=true, resistivity_model::NeoResistivityModel=SauterNeoModel(), @@ -247,17 +247,17 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; q, zeff; lnLamb=lnL) end - # dr_val: per-surface resistive interchange index D_R = E + F + H² + # delta_crit_D_R: per-surface resistive interchange index D_R = E + F + H² # (Glasser-Greene-Johnson 1975). Used by `_solve_dc_tmp` to compute # the χ_‖-matching critical-Δ via Connor-Hastie-Helander 2015 Eq. 59, # which has `(−D_R)` as a multiplier. NOT the Mercier index # D_I = E + F + H − 1/4 (see this function's docstring); we use the # physically correct D_R here. - dr_val_k = if dr_val === nothing + delta_crit_D_R_k = if delta_crit_D_R === nothing rg === nothing && throw( ArgumentError( - "build_slayer_inputs: dr_val=nothing " * + "build_slayer_inputs: delta_crit_D_R=nothing " * "requires `sing.restype` populated by " * "ForceFreeStates.resist_eval_all!. " * "Surface k=$k has restype=nothing." @@ -265,19 +265,19 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; ) rg.E + rg.F + rg.H^2 else - _eval(dr_val, psi) + _eval(delta_crit_D_R, psi) end - # dgeo_val: only used by dc_type=:toroidal (the Connor-Hastie- + # delta_crit_geo_factor: only used by delta_crit_type=:toroidal (the Connor-Hastie- # Helander 2015 formula). Auto-derivation requires ⟨|∇ψ|²⟩ FSA # which the current `ResistGeometry` doesn't expose; for now we - # require an explicit value if the toroidal dc_type is selected. - dgeo_val_k = if dgeo_val === nothing - dc_type === :toroidal && + # require an explicit value if the toroidal delta_crit_type is selected. + delta_crit_geo_factor_k = if delta_crit_geo_factor === nothing + delta_crit_type === :toroidal && throw( ArgumentError( - "build_slayer_inputs: dc_type=:toroidal " * - "needs `dgeo_val` (Connor 2015 PPCF 57 " * + "build_slayer_inputs: delta_crit_type=:toroidal " * + "needs `delta_crit_geo_factor` (Connor 2015 PPCF 57 " * "065001 Eq. 59 geometric factor). " * "Auto-derivation from equilibrium not " * "yet implemented; pass a scalar / vector " * @@ -286,7 +286,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; ) 0.0 else - _eval(dgeo_val, psi) + _eval(delta_crit_geo_factor, psi) end out[k] = slayer_parameters(; @@ -297,9 +297,9 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; chi_perp=_eval(chi_perp, psi), chi_tor=_eval(chi_tor, psi), m=m_res, n=n_res, - dr_val=dr_val_k, - dgeo_val=dgeo_val_k, - dc_type=dc_type, ising=k, + delta_crit_D_R=delta_crit_D_R_k, + delta_crit_geo_factor=delta_crit_geo_factor_k, + delta_crit_type=delta_crit_type, ising=k, resistivity_model=resistivity_model, f_trap=f_trap_kw, nu_e_star=nu_e_star_kw, diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index 57ee30763..d24de21e1 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -40,12 +40,12 @@ 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) | +| `delta_crit_D_R` | Radial width parameter at surface (input to dc_tmp) | +| `delta_crit_geo_factor` | Geometric Δ (Shafranov shift factor) | | `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 | -| `dc_type` | Selector for `dc_tmp` formula | +| `delta_crit_type` | Selector for `dc_tmp` formula | The complex normalized growth rate `Q = ω + iγ` is **not** stored here; it is passed as a separate argument to `solve_inner`. @@ -77,19 +77,19 @@ Base.@kwdef struct SLAYERParameters <: InnerLayerParameters R0::Float64 bt::Float64 sval_r::Float64 - dr_val::Float64 = 0.0 - dgeo_val::Float64 = 0.0 + delta_crit_D_R::Float64 = 0.0 + delta_crit_geo_factor::Float64 = 0.0 eta::Float64 d_beta::Float64 # Critical-Δ offset dc_tmp::Float64 = 0.0 - dc_type::Symbol = :none + delta_crit_type::Symbol = :none end -# Allowed dc_type values for the critical-Δ offset. `:none` is the default +# Allowed delta_crit_type values for the critical-Δ offset. `:none` is the default # `dc_tmp = 0` branch. -const ALLOWED_DC_TYPES = (:none, :lar, :rfitzp, :toroidal) +const ALLOWED_DELTA_CRIT_TYPES = (:none, :lar, :fitzpatrick, :toroidal) """ r_based_shear(rs, q, dq_dpsi, da_dpsi) -> Float64 @@ -118,14 +118,14 @@ end # Internal: solve the Wd self-consistency loop for the chi_parallel-based # critical Δ (Connor-Hastie-Helander 2015). Returns dc_tmp as a Float64. -function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real, +function _solve_dc_tmp(; delta_crit_type::Symbol, delta_crit_D_R::Real, delta_crit_geo_factor::Real, chi_perp::Real, t_e::Real, zeff::Real, tau_ee::Real, rs::Real, R0::Real, sval_r::Real, n_tor::Integer, max_iter::Integer=100, tol::Real=1e-10) - dc_type in ALLOWED_DC_TYPES || - throw(ArgumentError("SLAYERParameters: unknown dc_type=$dc_type. " * - "Allowed: $(ALLOWED_DC_TYPES)")) - (dc_type === :none || dr_val == 0.0) && return 0.0 + delta_crit_type in ALLOWED_DELTA_CRIT_TYPES || + throw(ArgumentError("SLAYERParameters: unknown delta_crit_type=$delta_crit_type. " * + "Allowed: $(ALLOWED_DELTA_CRIT_TYPES)")) + (delta_crit_type === :none || delta_crit_D_R == 0.0) && return 0.0 vte = sqrt(2.0 * t_e * E_CHG / M_E) chi_par_smfp = (1.581 * tau_ee * vte^2) / (1.0 + 0.2535 * zeff) @@ -150,15 +150,15 @@ function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real, chi_par_lmfp = (2.0 * R0 * vte) / (sqrt(π) * n_tor * sval_r * Wd) chi_par = (chi_par_smfp * chi_par_lmfp) / (chi_par_smfp + chi_par_lmfp) - if dc_type === :lar - return 0.5 * (-dr_val) * π^1.5 * + if delta_crit_type === :lar + return 0.5 * (-delta_crit_D_R) * π^1.5 * (chi_par / chi_perp)^0.25 * sqrt((n_tor * sval_r) / (R0 * rs)) - elseif dc_type === :rfitzp - return -(sqrt(2.0) * π^1.5 * dr_val) / Wd - elseif dc_type === :toroidal - return 0.5 * (-dr_val) * π^1.5 * - (chi_par / chi_perp)^0.25 * dgeo_val + elseif delta_crit_type === :fitzpatrick + return -(sqrt(2.0) * π^1.5 * delta_crit_D_R) / Wd + elseif delta_crit_type === :toroidal + return 0.5 * (-delta_crit_D_R) * π^1.5 * + (chi_par / chi_perp)^0.25 * delta_crit_geo_factor end return 0.0 end @@ -168,8 +168,8 @@ end qval, sval_r, bt, rs, R0, mu_i, zeff, chi_perp, chi_tor, m, n, - dr_val=0.0, dgeo_val=0.0, - dc_type=:none, ising=0, + delta_crit_D_R=0.0, delta_crit_geo_factor=0.0, + delta_crit_type=:none, ising=0, resistivity_model=SauterNeoModel(), f_trap=nothing, nu_e_star=nothing, R_major_eff=nothing, @@ -199,8 +199,8 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl - `zeff` -- effective charge - `chi_perp`, `chi_tor` -- perpendicular / toroidal heat diffusivity [m²/s] - `m`, `n` -- poloidal / toroidal mode numbers at the surface - - `dr_val`, `dgeo_val` -- inputs for the critical-Δ formula - - `dc_type` -- one of `:none`, `:lar`, `:rfitzp`, `:toroidal` + - `delta_crit_D_R`, `delta_crit_geo_factor` -- inputs for the critical-Δ formula + - `delta_crit_type` -- one of `:none`, `:lar`, `:fitzpatrick`, `:toroidal` - `ising` -- singular-surface index for traceability # Resistivity kwargs @@ -246,8 +246,8 @@ function slayer_parameters(; rs::Real, R0::Real, mu_i::Real, zeff::Real, chi_perp::Real, chi_tor::Real, m::Integer, n::Integer, - dr_val::Real=0.0, dgeo_val::Real=0.0, - dc_type::Symbol=:none, ising::Integer=0, + delta_crit_D_R::Real=0.0, delta_crit_geo_factor::Real=0.0, + delta_crit_type::Symbol=:none, ising::Integer=0, resistivity_model::NeoResistivityModel=SauterNeoModel(), f_trap::Union{Real,Nothing}=nothing, nu_e_star::Union{Real,Nothing}=nothing, @@ -344,7 +344,7 @@ function slayer_parameters(; delta_n = lu^(1.0 / 3.0) / rs # Critical-Δ offset from chi_parallel matching - dc_tmp = _solve_dc_tmp(; dc_type=dc_type, dr_val=dr_val, dgeo_val=dgeo_val, + dc_tmp = _solve_dc_tmp(; delta_crit_type=delta_crit_type, delta_crit_D_R=delta_crit_D_R, delta_crit_geo_factor=delta_crit_geo_factor, chi_perp=chi_perp, t_e=t_e, zeff=zeff, tau_ee=tau_ee, rs=rs, R0=R0, sval_r=sval_r, n_tor=n) @@ -356,8 +356,8 @@ function slayer_parameters(; Q_e=Q_e, Q_i=Q_i, iota_e=iota_e, tauk=tauk, tau_r=tau_r, delta_n=delta_n, rs=rs, R0=R0, bt=bt, sval_r=sval_r, - dr_val=dr_val, dgeo_val=dgeo_val, + delta_crit_D_R=delta_crit_D_R, delta_crit_geo_factor=delta_crit_geo_factor, eta=eta, d_beta=d_beta, - dc_tmp=dc_tmp, dc_type=dc_type + dc_tmp=dc_tmp, delta_crit_type=delta_crit_type ) end diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 4175a62b4..dbfe528dd 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -23,7 +23,7 @@ Build the ψ-quadrature node list `[x0, interior points strictly inside (x0, xou resonances); this function owns the ordering: sort, drop near-duplicates (closer than `PANEL_MERGE_ATOL`, e.g. a kinetic resonance coinciding with a rational), and drop points within `PANEL_MERGE_ATOL` of a bound to avoid degenerate panels. Paneling the integral at -these surfaces puts the resonant torque-density peaks (reg_spot/collisionally broadened, but +these surfaces puts the resonant torque-density peaks (regularization_width/collisionally broadened, but narrow in ψ) on Gauss-Kronrod interval endpoints, which the rule handles natively instead of hunting them by adaptive bisection. """ diff --git a/src/KineticForces/EnergyIntegration.jl b/src/KineticForces/EnergyIntegration.jl index 36cc1ecef..191eb4bfb 100644 --- a/src/KineticForces/EnergyIntegration.jl +++ b/src/KineticForces/EnergyIntegration.jl @@ -73,12 +73,12 @@ For CGL there is no resonance denominator: N_cgl = x^2.5 / (i·n). x25 = x * x * sqrt(x) # x^2.5 fx = if p.f0type == "maxwellian" ComplexF64((p.we + p.wn + p.wt * (x - 1.5)) * x25) - elseif p.f0type == "jkp" + elseif p.f0type == "park" ComplexF64((p.we + p.wn + p.wt * 2) * x25) elseif p.f0type == "cgl" complex(0.0, -x25 / p.n) # x^2.5 / (i*n) else - error("f0type must be maxwellian, jkp, or cgl") + error("f0type must be maxwellian, park, or cgl") end if p.qt fx *= (x - 2.5) @@ -101,12 +101,12 @@ never carries a CGL numerator (CGL has no pole). a = p.we + p.wn + p.wt * (x - 1.5) nn = ComplexF64(a * x25) dn = ComplexF64(p.wt * x25 + a * 2.5 * x15) # d/dx[(…)·x^2.5] - elseif p.f0type == "jkp" + elseif p.f0type == "park" a = p.we + p.wn + p.wt * 2 nn = ComplexF64(a * x25) dn = ComplexF64(a * 2.5 * x15) else - error("_energy_numerator_deriv supports maxwellian and jkp") + error("_energy_numerator_deriv supports maxwellian and park") end return p.qt ? dn * (x - 2.5) + nn : dn # d/dx[N·(x-2.5)] = N′·(x-2.5) + N end @@ -302,7 +302,7 @@ all collisionalities: the collisionless case (ν ≡ 0) is the exact ν→0 limi its real-axis pole resolved analytically (see `_integrate_energy_resonant`). Collision operator types (`nutype`): `"zero"`, `"small"`, `"krook"`, `"harmonic"`. -Distribution function types (`f0type`): `"maxwellian"`, `"jkp"`, `"cgl"`. +Distribution function types (`f0type`): `"maxwellian"`, `"park"`, `"cgl"`. `ximag` is accepted for backward compatibility but no longer used — resonance poles are now handled analytically rather than by contour deformation. diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 6df40e6a8..c5a85e33a 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -126,7 +126,7 @@ ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForce # Energy integration parameters nutype::String = "harmonic" # Collision operator: "zero", "small", "krook", "harmonic" - f0type::String = "maxwellian" # Distribution function: "maxwellian", "jkp", "cgl" + f0type::String = "maxwellian" # Distribution function: "maxwellian", "park", "cgl" # Diagnostic parameters psilims::Vector{Float64} = [0.0, 1.0] # Integration limits in psi diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl index e491477c4..16d9d57a1 100644 --- a/src/PerturbedEquilibrium/FieldReconstruction.jl +++ b/src/PerturbedEquilibrium/FieldReconstruction.jl @@ -18,7 +18,7 @@ where χ₁ = 2π·Ψ₀ [Park Phys. Plasmas 14, 052110 (2007) eq. 8-10]. Clebsch displacement components for PENTRC (matches Fortran gpout_xclebsch): ξ^ψ = xsp_mn (unregularized) - ∂ξ^ψ/∂ψ = xmp1_mn (regularized: xsp1 * singfac²/(singfac² + reg_spot²)) + ∂ξ^ψ/∂ψ = xmp1_mn (regularized: xsp1 * singfac²/(singfac² + regularization_width²)) ξ^α = xms_mn (regularized: -A⁻¹(B·xmp1 + C·xsp), divided by χ₁ in output) Contravariant displacement from Jacobian convolution (matches Fortran gpeq_contra): @@ -56,14 +56,14 @@ Tuple of (xi_modes, b_modes) NamedTuples: - `xi_modes.clebsch_psi`: ξ^ψ for PENTRC (= xi_psi, unregularized) - `xi_modes.clebsch_psi1`: ∂ξ^ψ/∂ψ regularized for PENTRC - `xi_modes.clebsch_alpha`: ξ^α/χ₁ regularized for PENTRC (divided by χ₁ per gpout_xclebsch) - - `xi_modes.theta_reg`: ξ^θ regularized (= xmt, from gpeq_contra with reg_spot smoothing) - - `xi_modes.zeta_reg`: ξ^ζ regularized (= xmz, from gpeq_contra with reg_spot smoothing) + - `xi_modes.theta_reg`: ξ^θ regularized (= xmt, from gpeq_contra with regularization_width smoothing) + - `xi_modes.zeta_reg`: ξ^ζ regularized (= xmz, from gpeq_contra with regularization_width smoothing) - `xi_modes.cova_psi/theta/zeta`: covariant displacement (from gpeq_cova) - `b_modes.psi`: b^ψ [npsi, mpert] - `b_modes.b_psi_area_weighted`: b^ψ / ⟨J·|∇ψ|⟩_θ (area-normalized, for b_n computation) - `b_modes.theta`: b^θ [npsi, mpert] - `b_modes.zeta`: b^ζ [npsi, mpert] - - `b_modes.theta_reg/zeta_reg`: regularized b^θ, b^ζ (from gpeq_sol with reg_spot smoothing) + - `b_modes.theta_reg/zeta_reg`: regularized b^θ, b^ζ (from gpeq_sol with regularization_width smoothing) - `b_modes.cova_psi/theta/zeta`: covariant field (from gpeq_cova) """ function reconstruct_physical_fields( @@ -177,8 +177,8 @@ function reconstruct_physical_fields( psi_J=xwp_modes, # J·ξ^ψ (Jacobian-weighted, from gpeq_contra) theta=xwt_modes, # ξ^θ contravariant (from gpeq_contra) zeta=xwz_modes, # ξ^ζ contravariant (from gpeq_contra) - theta_reg=xmt_modes, # ξ^θ regularized (from gpeq_contra, smoothed by reg_spot) - zeta_reg=xmz_modes, # ξ^ζ regularized (from gpeq_contra, smoothed by reg_spot) + theta_reg=xmt_modes, # ξ^θ regularized (from gpeq_contra, smoothed by regularization_width) + zeta_reg=xmz_modes, # ξ^ζ regularized (from gpeq_contra, smoothed by regularization_width) clebsch_psi=clebsch_psi, # ξ^ψ for PENTRC clebsch_psi1=clebsch_psi1, # ∂ξ^ψ/∂ψ regularized for PENTRC clebsch_alpha=clebsch_alpha, # ξ^α/χ₁ regularized for PENTRC @@ -192,7 +192,7 @@ function reconstruct_physical_fields( b_psi_area_weighted=Jb_psi_modes, # b^ψ / ⟨J·|∇ψ|⟩_θ (area-normalized, for b_n) theta=b_theta_modes, # b^θ unregularized zeta=b_zeta_modes, # b^ζ unregularized - theta_reg=b_theta_reg, # b^θ regularized (from gpeq_sol with reg_spot) + theta_reg=b_theta_reg, # b^θ regularized (from gpeq_sol with regularization_width) zeta_reg=b_zeta_reg, # b^ζ regularized cova_psi=bvp_modes, # covariant b_ψ (from gpeq_cova) cova_theta=bvt_modes, # covariant b_θ (from gpeq_cova) @@ -329,10 +329,10 @@ Compute Clebsch displacement components for PENTRC output. Matches Fortran gpeq_sol regularization + gpout_xclebsch output convention: - `clebsch_psi` = ξ^ψ (unregularized, same as xi_psi_modes) - - `clebsch_psi1` = xmp1 = ∂ξ^ψ/∂ψ × singfac²/(singfac² + reg_spot²) + - `clebsch_psi1` = xmp1 = ∂ξ^ψ/∂ψ × singfac²/(singfac² + regularization_width²) - `clebsch_alpha` = xms/χ₁ (regularized ξ^α divided by χ₁ per gpout_xclebsch convention) -When reg_spot=0, clebsch_psi1 = xi_psi1 and clebsch_alpha = xi_s/χ₁ (no regularization). +When regularization_width=0, clebsch_psi1 = xi_psi1 and clebsch_alpha = xi_s/χ₁ (no regularization). The regularized xms is computed as -A⁻¹(B·xmp1 + C·xsp) matching Fortran gpeq_sol, where A, B, C are the stability matrices evaluated at each ψ via ffit interpolants. @@ -357,10 +357,10 @@ function compute_clebsch_displacements( clebsch_psi1 = copy(xi_psi1_modes) # will be regularized below clebsch_alpha = xi_s_modes ./ chi1 # ξ^α/χ₁ (will be regularized below) - reg_spot = ctrl.reg_spot - @assert reg_spot >= 0 "reg_spot must be non-negative (got $reg_spot)" + regularization_width = ctrl.regularization_width + @assert regularization_width >= 0 "regularization_width must be non-negative (got $regularization_width)" - if reg_spot == 0 + if regularization_width == 0 return clebsch_psi, clebsch_psi1, clebsch_alpha end @@ -390,7 +390,7 @@ function compute_clebsch_displacements( for ipert in 1:mpert m = mlow + ipert - 1 singfac = m - nn * q - reg_factor = singfac^2 / (singfac^2 + reg_spot^2) + reg_factor = singfac^2 / (singfac^2 + regularization_width^2) clebsch_psi1[ipsi, ipert] = xi_psi1_modes[ipsi, ipert] * reg_factor xmp1_vec[ipert] = clebsch_psi1[ipsi, ipert] end @@ -498,7 +498,7 @@ function compute_contra_displacements( mlow = ffs_intr.mlow nn = ffs_intr.nlow chi1 = 2π * equil.psio - reg_spot = ctrl.reg_spot + regularization_width = ctrl.regularization_width fc = metric.fourier_coeffs xwp_modes = zeros(ComplexF64, npsi, mpert) @@ -581,14 +581,14 @@ function compute_contra_displacements( # Regularize xwt/xwz → xmt/xmz (matches Fortran gpeq_contra) xmt_modes = copy(xwt_modes) xmz_modes = copy(xwz_modes) - if reg_spot > 0 + if regularization_width > 0 Threads.@threads :static for ipsi in 1:npsi psi_norm = psi_grid[ipsi] q = equil.profiles.q_spline(psi_norm) for ipert in 1:mpert m = mlow + ipert - 1 singfac = m - nn * q - reg_factor = singfac^2 / (singfac^2 + reg_spot^2) + reg_factor = singfac^2 / (singfac^2 + regularization_width^2) xmt_modes[ipsi, ipert] = xwt_modes[ipsi, ipert] * reg_factor xmz_modes[ipsi, ipert] = xwz_modes[ipsi, ipert] * reg_factor end diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index 5e9429d5b..3b83889c2 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -116,7 +116,7 @@ function compute_perturbed_equilibrium( # Same control surface as the coil branch above: psilim, the integration # limit (Fortran: gpec/gpec.f:431 `field_bs_psi(psilim, ...)`). Without it # the normalization was taken on the equilibrium-spline limit, which differs - # whenever dmlim/qhigh/psiedge truncation moves psilim inward. + # whenever dmlim/qhigh/dW_edge_scan_start truncation moves psilim inward. convert_forcing_normalization!(modes_n, norm_tag, equil, n, minimum(m_vals), maximum(m_vals); psi=ffs_intr.psilim) end diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index 96f982386..775b46b90 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -27,7 +27,7 @@ Medium Priority (defer for MWE): Regularization: # High Priority (MWE) - - `reg_spot::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0. + - `regularization_width::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0. """ @kwdef mutable struct PerturbedEquilibriumControl # High Priority (MWE) @@ -47,7 +47,7 @@ Regularization: # Regularization width for singular surface smoothing (matches Fortran gpec.f reg_spot). # Set to 0 to disable regularization. Must be non-negative. - reg_spot::Float64 = 5e-2 + regularization_width::Float64 = 5e-2 end """ diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl index 322c6cf9c..f563d4c2a 100644 --- a/src/PerturbedEquilibrium/Utils.jl +++ b/src/PerturbedEquilibrium/Utils.jl @@ -260,7 +260,7 @@ const PE_H5_ANNOTATIONS = [ "Response/xi_cova_theta" => (; long_name="covariant poloidal displacement ξ_θ", units="m^2", dims=("psi", "mode")), "Response/xi_cova_zeta" => (; long_name="covariant toroidal displacement ξ_ζ", units="m^2", dims=("psi", "mode")), "Response/clebsch_psi" => (; long_name="Clebsch displacement component ξ^ψ (PENTRC input, gpout_xclebsch convention)", dims=("psi", "mode")), - "Response/clebsch_psi1" => (; long_name="regularized ψ_N derivative of ξ^ψ (× singfac²/(singfac²+reg_spot²))", dims=("psi", "mode")), + "Response/clebsch_psi1" => (; long_name="regularized ψ_N derivative of ξ^ψ (× singfac²/(singfac²+regularization_width²))", dims=("psi", "mode")), "Response/clebsch_alpha" => (; long_name="Clebsch displacement component ξ^α/χ₁ (PENTRC input, gpout_xclebsch convention)", dims=("psi", "mode")), "Response/xi_n" => (; long_name="physical normal displacement ξ_n", units="m", dims=("psi", "mode")), "Response/xi_R" => (; long_name="cylindrical displacement component ξ_R (mode space)", units="m", dims=("psi", "mode")), diff --git a/src/Rerun.jl b/src/Rerun.jl index 4eab2a985..9455c8001 100644 --- a/src/Rerun.jl +++ b/src/Rerun.jl @@ -265,6 +265,7 @@ function build_inputs_from_h5(args::Vector{String}) " source: $(abspath(source_h5))\n" * " output: $(abspath(joinpath(output_dir, output_name)))\n$_BANNER" + _rename_keys!(inputs["Equilibrium"], _RENAMED_EQUIL_KEYS, "Equilibrium") _drop_deprecated_keys!(inputs["Equilibrium"], _DEPRECATED_EQUIL_KEYS, "Equilibrium") eq_config = Equilibrium.EquilibriumConfig(inputs["Equilibrium"], output_dir) # Clear eq_filename: unused on replay, and a stale absolute path could mislead downstream code. diff --git a/src/Tearing/Dispersion/CoupledFullMatch.jl b/src/Tearing/Dispersion/CoupledFullMatch.jl index 8bf5fb97c..6360c8fbe 100644 --- a/src/Tearing/Dispersion/CoupledFullMatch.jl +++ b/src/Tearing/Dispersion/CoupledFullMatch.jl @@ -94,7 +94,7 @@ end Construct the 4m × 4m dispersion matrix driver. `dp_raw` must be the 2m × 2m matrix in side-major ordering (the `intr.delta_prime_raw` field populated by `ForceFreeStates.compute_delta_prime_matrix!` on the -`use_parallel=true` path). `rotation[k]` is the per-surface rotation +`integrator="stride"` path). `rotation[k]` is the per-surface rotation frequency; it shifts the per-surface inner Q argument by `i·ntor·rotation[k]`. Default zero rotation matches the static-equilibrium case. diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 17128a3d6..49b5518b0 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -21,8 +21,8 @@ constructor. - `scan_mode` -- `:amr` (default) or `:brute_force` - `coupling_mode` -- `:uncoupled` (default, per-surface) or `:coupled` (multi-surface determinant) - - `dc_type` -- critical-Δ offset selector, one of `:none`, `:lar`, - `:rfitzp`, `:toroidal` (χ_‖-matching critical-Δ formulas, + - `delta_crit_type` -- critical-Δ offset selector, one of `:none`, `:lar`, + `:fitzpatrick`, `:toroidal` (χ_‖-matching critical-Δ formulas, Connor-Hastie-Helander 2015) - `msing_max` -- number of surfaces to include in the coupled determinant (default 3; capped at `length(sings)` at runtime) @@ -36,10 +36,10 @@ constructor. diffusivity [m²/s], used only when the kinetic file carries no usable `chi_e`/`chi_phi` profile (dataset absent or all-zero); otherwise the file's χ⊥(ψ)/χ_φ(ψ) take precedence - - `dr_val`, `dgeo_val` -- critical-Δ formula inputs. `nothing` (default) - auto-derives them from the equilibrium: `dr_val` from the resistive - interchange index `D_R = E + F + H²` at each surface, `dgeo_val` from the - toroidal geometric factor (required only by `dc_type=:toroidal`). Supply a + - `delta_crit_D_R`, `delta_crit_geo_factor` -- critical-Δ formula inputs. `nothing` (default) + auto-derives them from the equilibrium: `delta_crit_D_R` from the resistive + interchange index `D_R = E + F + H²` at each surface, `delta_crit_geo_factor` from the + toroidal geometric factor (required only by `delta_crit_type=:toroidal`). Supply a scalar only to override the auto-derivation; an explicit `0.0` disables the critical-Δ offset (Δ_crit ≡ 0) - `theta_sample` -- poloidal angle at which to sample minor radius @@ -99,7 +99,7 @@ there is one consistent interface for resistive and kinetic profiles. inner_model::Symbol = :slayer_fitzpatrick scan_mode::Symbol = :amr coupling_mode::Symbol = :uncoupled - dc_type::Symbol = :none + delta_crit_type::Symbol = :none msing_max::Int = 3 bt::Union{Float64,Nothing} = nothing @@ -107,8 +107,8 @@ there is one consistent interface for resistive and kinetic profiles. zeff::Float64 = 1.0 chi_perp::Float64 = 1.0 chi_tor::Float64 = 1.0 - dr_val::Union{Float64,Nothing} = nothing - dgeo_val::Union{Float64,Nothing} = nothing + delta_crit_D_R::Union{Float64,Nothing} = nothing + delta_crit_geo_factor::Union{Float64,Nothing} = nothing theta_sample::Float64 = 0.0 resistivity_model::Symbol = :sauter lnLambda_form::Symbol = :nrl @@ -158,7 +158,7 @@ end const _VALID_INNER_MODELS = (:slayer_fitzpatrick, :ggj_shooting, :ggj_galerkin) const _VALID_SCAN_MODES = (:amr, :brute_force) const _VALID_COUPLING_MODES = (:uncoupled, :coupled) -const _VALID_DC_TYPES = (:none, :lar, :rfitzp, :toroidal) +const _VALID_DELTA_CRIT_TYPES = (:none, :lar, :fitzpatrick, :toroidal) const _VALID_RESISTIVITY_MODELS = (:sauter, :redl, :spitzer, :spitzer_harm) const _VALID_LNLAMBDA_FORMS = (:nrl, :sauter, :wesson) @@ -172,9 +172,9 @@ function validate(ctrl::SLAYERControl) ctrl.coupling_mode in _VALID_COUPLING_MODES || throw(ArgumentError("SLAYERControl: coupling_mode=$(ctrl.coupling_mode) " * "not in $(_VALID_COUPLING_MODES)")) - ctrl.dc_type in _VALID_DC_TYPES || - throw(ArgumentError("SLAYERControl: dc_type=$(ctrl.dc_type) " * - "not in $(_VALID_DC_TYPES)")) + ctrl.delta_crit_type in _VALID_DELTA_CRIT_TYPES || + throw(ArgumentError("SLAYERControl: delta_crit_type=$(ctrl.delta_crit_type) " * + "not in $(_VALID_DELTA_CRIT_TYPES)")) ctrl.resistivity_model in _VALID_RESISTIVITY_MODELS || throw(ArgumentError("SLAYERControl: resistivity_model=$(ctrl.resistivity_model) " * "not in $(_VALID_RESISTIVITY_MODELS)")) @@ -227,6 +227,20 @@ function slayer_control_from_toml(section::AbstractDict) end end + # Deprecated [SLAYER] spellings: renamed keys and the :rfitzp value alias, accepted + # with a warning until removal after v2.0.0. + for (old, new) in ("dc_type" => "delta_crit_type", "dr_val" => "delta_crit_D_R", + "dgeo_val" => "delta_crit_geo_factor") + haskey(flat, old) || continue + @warn "`$old` in [SLAYER] was renamed to `$new`; the old key is deprecated and will be removed after v2.0.0." + haskey(flat, new) || (flat[new] = flat[old]) + delete!(flat, old) + end + if get(flat, "delta_crit_type", "") in ("rfitzp", :rfitzp) + @warn "`delta_crit_type = \"rfitzp\"` in [SLAYER] is deprecated; use `delta_crit_type = \"fitzpatrick\"`. The old value will be removed after v2.0.0." + flat["delta_crit_type"] = "fitzpatrick" + end + # Validate keys against the struct fields field_names = Set(String.(fieldnames(SLAYERControl))) unknown = [k for k in keys(flat) if !(k in field_names)] @@ -239,12 +253,12 @@ function slayer_control_from_toml(section::AbstractDict) kwargs = Dict{Symbol,Any}() for (k, v) in flat sym = Symbol(k) - if sym in (:inner_model, :scan_mode, :coupling_mode, :dc_type, + if sym in (:inner_model, :scan_mode, :coupling_mode, :delta_crit_type, :resistivity_model, :lnLambda_form) kwargs[sym] = v isa Symbol ? v : Symbol(String(v)) elseif sym in (:Q_re_range, :Q_im_range) kwargs[sym] = _as_range(v) - elseif sym in (:bt, :dr_val, :dgeo_val) + elseif sym in (:bt, :delta_crit_D_R, :delta_crit_geo_factor) # Allow explicit nothing (auto-derive) or a number (override) kwargs[sym] = v === nothing ? nothing : Float64(v) elseif sym === :boxes diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index efff4ef0d..d081d49ee 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -162,12 +162,15 @@ function _write_per_surface!(g, params::AbstractVector{SLAYERParameters}, 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, + :rs, :R0, :bt, :sval_r, :eta, :d_beta, :dc_tmp) ps[String(fname)] = Float64[getfield(p, fname) for p in params] end - # Store dc_type per-surface as string array - ps["dc_type"] = String[String(p.dc_type) for p in params] + # HDF5 leaf names keep the legacy spellings (schema stability); the struct fields were renamed. + ps["dr_val"] = Float64[getfield(p, :delta_crit_D_R) for p in params] + ps["dgeo_val"] = Float64[getfield(p, :delta_crit_geo_factor) for p in params] + # Store the per-surface critical-Δ prescription label as string array + ps["dc_type"] = String[String(p.delta_crit_type) for p in params] # Full Δ' matrix, split real/imag dp = create_group(ps, "DpMatrix") diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 0e254e625..ea82968aa 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -396,9 +396,9 @@ function run_slayer(equil, ffs_intr, control::SLAYERControl; zeff=control.zeff, chi_perp=chi_perp, chi_tor=chi_tor, - dr_val=control.dr_val, - dgeo_val=control.dgeo_val, - dc_type=control.dc_type, + delta_crit_D_R=control.delta_crit_D_R, + delta_crit_geo_factor=control.delta_crit_geo_factor, + delta_crit_type=control.delta_crit_type, theta=control.theta_sample, resistivity_model=_build_resistivity_model(control.resistivity_model), lnLambda_form=control.lnLambda_form) diff --git a/test/runtests.jl b/test/runtests.jl index 76eca69d1..7368d3b38 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -44,6 +44,7 @@ else include("./runtests_dispersion_amr.jl") include("./runtests_dispersion_polish.jl") include("./runtests_slayer_runner.jl") + include("./runtests_toml_backcompat.jl") include("./runtests_kinetic.jl") include("./runtests_fullruns.jl") include("./runtests_coils.jl") diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index dd1e49ce5..5fb0e22ac 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -10,7 +10,7 @@ eq_filename=joinpath(data_dir, "EQDSK_COCOS_02"), eq_type="efit", jac_type="boozer", - grid_type="ldp", + grid_type="rational_packed", psilow=0.01, psihigh=0.994 ) @@ -25,7 +25,7 @@ eq_filename=joinpath(data_dir, "EQDSK_COCOS_02"), eq_type="efit_arclength", jac_type="boozer", - grid_type="ldp", + grid_type="rational_packed", psilow=0.01, psihigh=0.994 ) @@ -49,7 +49,7 @@ eq_filename=joinpath(data_dir, "EQDSK_COCOS_02"), eq_type="efit_by_inversion", jac_type="boozer", - grid_type="ldp", + grid_type="rational_packed", psilow=0.01, psihigh=0.994 ) @@ -90,7 +90,7 @@ eq_filename=joinpath(data_dir, "INP1_binary"), eq_type="chease_binary", jac_type="boozer", - grid_type="ldp", + grid_type="rational_packed", psilow=0.01, psihigh=0.994, r0exp=6.8, @@ -106,7 +106,7 @@ eq_filename=joinpath(data_dir, "INP1_ascii"), eq_type="chease_ascii", jac_type="boozer", - grid_type="ldp", + grid_type="rational_packed", psilow=0.01, psihigh=0.994, r0exp=6.8, @@ -201,7 +201,7 @@ lar_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(; eq_type="lar", jac_type="boozer", - grid_type="ldp", + grid_type="rational_packed", psilow=0.01, psihigh=0.99 ) @@ -367,7 +367,7 @@ mpsi=64, mtheta=128) eq_config = Eq.EquilibriumConfig(; eq_type="sol", eq_filename="unused", - jac_type="pest", grid_type="ldp", + jac_type="pest", grid_type="rational_packed", psilow=1e-4, psihigh=0.99999, mpsi=mpsi, mtheta=mtheta) sol_config = Eq.SolovevConfig(64, 64, 64, e, a, r0, q0, 1.0, 1.0, 1.0) dri = Eq.sol_run(eq_config, sol_config) diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 7adb75203..27e6804c3 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -432,7 +432,7 @@ end example_dir = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") inputs = TOML.parsefile(joinpath(example_dir, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_parallel"] = false + inputs["ForceFreeStates"]["integrator"] = "serial" inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false intr = FFS.ForceFreeStatesInternal(; dir_path=example_dir) ctrl = FFS.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index d23a00790..4e672a63b 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -227,14 +227,14 @@ using TOML # The energy eigenvalue et[1] should match to within 2%. # # Bidirectional FM integration (crossing chunks integrated backward) is the - # default for use_parallel=true. It keeps FM propagators well-conditioned for + # default for integrator="stride". It keeps FM propagators well-conditioned for # both small-N (Solovev N=8, tested here) and large-N (DIIID N=26, tested below). ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") - function run_solovev(use_parallel) + function run_solovev(integrator) inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_parallel"] = use_parallel + inputs["ForceFreeStates"]["integrator"] = integrator intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=ex) ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) @@ -258,8 +258,8 @@ using TOML return real(vac.et[1]), intr end - et_std, intr_std = run_solovev(false) - et_par, intr_par = run_solovev(true) + et_std, intr_std = run_solovev("serial") + et_par, intr_par = run_solovev("stride") # Energy eigenvalue matches to 2% @test isapprox(et_par, et_std; rtol=0.02) @@ -279,10 +279,10 @@ using TOML # This is the key regression test for the bidirectional parallel FM fix. ex = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example") - function run_diiid(use_parallel) + function run_diiid(integrator) inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_parallel"] = use_parallel + inputs["ForceFreeStates"]["integrator"] = integrator inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=ex) ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; @@ -317,7 +317,7 @@ using TOML return real(vac.et[1]), intr end - et_par, intr_par = run_diiid(true) + et_par, intr_par = run_diiid("stride") # Parallel FM et[1] regression — pinned tightly, NOT bracketed. et[1] is grid- and # equilibrium-sensitive (auto-mpsi gives a spurious value; a wrong grid/Ip shifts it), so @@ -386,12 +386,12 @@ using TOML # physically meaningful. BVP Δ' regression is concentrated on the DIIID-like # fixture below (intrinsically stable, well-conditioned BVP Δ'). - @testset "ξ functions bit-identical between use_parallel modes (populate_dense_xi)" begin - # When `ctrl.use_parallel = true` and `ctrl.populate_dense_xi = true` + @testset "ξ functions bit-identical between integrator modes (populate_dense_xi)" begin + # When `ctrl.integrator = "stride"` and `ctrl.populate_dense_xi = true` # (default), `parallel_eulerlagrange_integration` appends a serial # Euler-Lagrange pass and returns that fresh `odet` instead of the # propagator-BVP one. That dense pass invokes the SAME - # `eulerlagrange_integration` code path the serial `use_parallel = false` + # `eulerlagrange_integration` code path the serial `integrator = "serial"` # benchmark goes through with the SAME `(ctrl, equil, ffit, intr)` # inputs (BVP-only state on `intr` saved/restored across the pass), so # the resulting `psi_store` / `q_store` / `u_store` / `du_store` / @@ -403,10 +403,10 @@ using TOML # Run on both the small-N Solovev case and the large-N DIIID-like case # to catch any (m, IC, ψ)-dependent regression. - function run_and_capture(example_dir, use_parallel; populate_dense_xi=true) + function run_and_capture(example_dir, integrator; populate_dense_xi=true) inputs = TOML.parsefile(joinpath(example_dir, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_parallel"] = use_parallel + inputs["ForceFreeStates"]["integrator"] = integrator inputs["ForceFreeStates"]["populate_dense_xi"] = populate_dense_xi inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=example_dir) @@ -457,15 +457,15 @@ using TOML @testset "Solovev (small N)" begin ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") - odet_std = run_and_capture(ex, false) - odet_par = run_and_capture(ex, true; populate_dense_xi=true) + odet_std = run_and_capture(ex, "serial") + odet_par = run_and_capture(ex, "stride"; populate_dense_xi=true) assert_bit_identical(odet_std, odet_par) end @testset "DIIID-like (large N)" begin ex = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example") - odet_std = run_and_capture(ex, false) - odet_par = run_and_capture(ex, true; populate_dense_xi=true) + odet_std = run_and_capture(ex, "serial") + odet_par = run_and_capture(ex, "stride"; populate_dense_xi=true) assert_bit_identical(odet_std, odet_par) end @@ -477,8 +477,8 @@ using TOML # test above is meaningful — it's NOT trivially passing because # both modes accidentally produce the same sparse data. ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") - odet_std = run_and_capture(ex, false) - odet_sparse = run_and_capture(ex, true; populate_dense_xi=false) + odet_std = run_and_capture(ex, "serial") + odet_sparse = run_and_capture(ex, "stride"; populate_dense_xi=false) @test odet_sparse.step < odet_std.step @test length(odet_sparse.psi_store) < length(odet_std.psi_store) # The sparse solution is in the Riccati basis, so the derivative stores cannot be @@ -498,7 +498,7 @@ using TOML ex = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example") inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_parallel"] = true + inputs["ForceFreeStates"]["integrator"] = "stride" inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=ex) ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; diff --git a/test/runtests_rerun_from_h5.jl b/test/runtests_rerun_from_h5.jl index 2e364f2df..67c816d07 100644 --- a/test/runtests_rerun_from_h5.jl +++ b/test/runtests_rerun_from_h5.jl @@ -250,7 +250,7 @@ end Equil = GeneralizedPerturbedEquilibrium.Equilibrium config = Equil.EquilibriumConfig(; eq_filename=joinpath(@__DIR__, "test_data", "CHEASE_test_data", "INP1_ascii"), - eq_type="chease_ascii", jac_type="boozer", grid_type="ldp", + eq_type="chease_ascii", jac_type="boozer", grid_type="rational_packed", psilow=0.01, psihigh=0.994, r0exp=6.8, b0exp=7.4 ) src = Equil.read_chease_ascii(config) diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index 5e2c5a34f..7fb0925e1 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -67,13 +67,13 @@ @testset "build_slayer_inputs: returns correct per-surface data" begin sings = [_mk_sing(psi=0.3, q=2.0, q1=1.5, m=2, n=1), _mk_sing(psi=0.6, q=3.0, q1=2.5, m=3, n=1)] - # dr_val=0.0 bypasses the build_slayer_inputs requirement that sing.restype be + # delta_crit_D_R=0.0 bypasses the build_slayer_inputs requirement that sing.restype be # pre-populated by ForceFreeStates.resist_eval_all! — the test sings here are - # minimal stubs without restype, so we supply dr_val explicitly. + # minimal stubs without restype, so we supply delta_crit_D_R explicitly. # compute_omega_star=false makes Q_e/Q_i pass through directly from profiles.omega_e/i # rather than being recomputed from n_e/T_e/T_i gradients — required for the Q_e == # -tauk·omega_e(ψ) identity check below. - sl = build_slayer_inputs(equil, sings, profiles; bt=2.0, dr_val=0.0, + sl = build_slayer_inputs(equil, sings, profiles; bt=2.0, delta_crit_D_R=0.0, compute_omega_star=false) @test length(sl) == 2 @@ -110,43 +110,43 @@ @testset "build_slayer_inputs: chi_perp/chi_tor as scalars and callables" begin sings = [_mk_sing(psi=0.5, q=2.4, q1=1.2, m=2, n=1)] - # Scalar (dr_val=0.0 bypasses the sing.restype requirement; see comment above) + # Scalar (delta_crit_D_R=0.0 bypasses the sing.restype requirement; see comment above) sl_s = build_slayer_inputs(equil, sings, profiles; - bt=2.0, chi_perp=2.0, chi_tor=1.5, dr_val=0.0) + bt=2.0, chi_perp=2.0, chi_tor=1.5, delta_crit_D_R=0.0) # Callable with matching value chi_p(psi) = 2.0 + 0.0*psi chi_t(psi) = 1.5 + 0.0*psi sl_c = build_slayer_inputs(equil, sings, profiles; - bt=2.0, chi_perp=chi_p, chi_tor=chi_t, dr_val=0.0) + bt=2.0, chi_perp=chi_p, chi_tor=chi_t, delta_crit_D_R=0.0) @test sl_s[1].P_perp ≈ sl_c[1].P_perp @test sl_s[1].P_tor ≈ sl_c[1].P_tor # Callable with ψ-dependence changes the result chi_p_var(psi) = 1.0 + 10.0 * psi # χ⊥(0.5) = 6.0 > 2.0 sl_var = build_slayer_inputs(equil, sings, profiles; - bt=2.0, chi_perp=chi_p_var, chi_tor=1.5, dr_val=0.0) + bt=2.0, chi_perp=chi_p_var, chi_tor=1.5, delta_crit_D_R=0.0) # P_perp = τ_r · χ⊥ / r² grows with χ⊥, so the varying-χ case at # ψ=0.5 (χ⊥=6) gives a *larger* P_perp than the scalar χ⊥=2. @test sl_var[1].P_perp > sl_s[1].P_perp @test sl_var[1].P_perp ≈ sl_s[1].P_perp * 6.0 / 2.0 rtol = 1e-10 end - @testset "build_slayer_inputs: dc_type propagates and dr_val activates offset" begin + @testset "build_slayer_inputs: delta_crit_type propagates and delta_crit_D_R activates offset" begin sings = [_mk_sing(psi=0.5, q=2.4, q1=1.2, m=2, n=1)] - # dc_type=:none and dr_val=0.0 → dc_tmp = 0 regardless of dr_val + # delta_crit_type=:none and delta_crit_D_R=0.0 → dc_tmp = 0 regardless of delta_crit_D_R sl_none = build_slayer_inputs(equil, sings, profiles; - bt=2.0, dc_type=:none, dr_val=0.0) + bt=2.0, delta_crit_type=:none, delta_crit_D_R=0.0) @test sl_none[1].dc_tmp == 0.0 - # dc_type=:rfitzp with dr_val = 0 still gives zero + # delta_crit_type=:fitzpatrick with delta_crit_D_R = 0 still gives zero sl_rf0 = build_slayer_inputs(equil, sings, profiles; - bt=2.0, dc_type=:rfitzp, dr_val=0.0) + bt=2.0, delta_crit_type=:fitzpatrick, delta_crit_D_R=0.0) @test sl_rf0[1].dc_tmp == 0.0 - # dc_type=:rfitzp with dr_val > 0 → nonzero negative offset + # delta_crit_type=:fitzpatrick with delta_crit_D_R > 0 → nonzero negative offset sl_rf = build_slayer_inputs(equil, sings, profiles; - bt=2.0, dc_type=:rfitzp, dr_val=0.01) + bt=2.0, delta_crit_type=:fitzpatrick, delta_crit_D_R=0.01) @test sl_rf[1].dc_tmp < 0 @test isfinite(sl_rf[1].dc_tmp) end diff --git a/test/runtests_slayer_params.jl b/test/runtests_slayer_params.jl index 330ba7297..d1715bd38 100644 --- a/test/runtests_slayer_params.jl +++ b/test/runtests_slayer_params.jl @@ -6,7 +6,7 @@ # Reference inputs: a simple deuterium plasma case suitable for # hand-checking the SLAYER params formulas. - function _ref_kwargs(; dr_val=0.0, dc_type=:none) + function _ref_kwargs(; delta_crit_D_R=0.0, delta_crit_type=:none) return ( n_e=5.0e19, t_e=1000.0, t_i=1000.0, omega=0.0, omega_e=1.0e4, omega_i=5.0e3, @@ -14,7 +14,7 @@ rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0, chi_perp=1.0, chi_tor=1.0, m=2, n=1, - dr_val=dr_val, dgeo_val=0.5, dc_type=dc_type, + delta_crit_D_R=delta_crit_D_R, delta_crit_geo_factor=0.5, delta_crit_type=delta_crit_type, ising=3 ) end @@ -31,8 +31,8 @@ @test p.R0 == 1.7 @test p.bt == 2.0 @test p.sval_r == 1.0 - @test p.dc_tmp == 0.0 # dr_val == 0 ⇒ no offset - @test p.dc_type === :none + @test p.dc_tmp == 0.0 # delta_crit_D_R == 0 ⇒ no offset + @test p.delta_crit_type === :none # Trivially exact ratios @test p.tau ≈ 1.0 @@ -91,34 +91,34 @@ @test p.delta_n ≈ p.lu^(1 / 3) / p.rs rtol = 1e-12 end - @testset "Test 1b: dc_tmp formulas activate when dr_val ≠ 0" begin - # All four dc_type branches must produce finite, non-NaN values + @testset "Test 1b: dc_tmp formulas activate when delta_crit_D_R ≠ 0" begin + # All four delta_crit_type branches must produce finite, non-NaN values # and respect the signs/structure of the formulas in # the SLAYER params dc_tmp formulas. - p_none = slayer_parameters(; _ref_kwargs(; dr_val=0.01, dc_type=:none)...) - @test p_none.dc_tmp == 0.0 # :none ignores dr_val + p_none = slayer_parameters(; _ref_kwargs(; delta_crit_D_R=0.01, delta_crit_type=:none)...) + @test p_none.dc_tmp == 0.0 # :none ignores delta_crit_D_R - p_lar = slayer_parameters(; _ref_kwargs(; dr_val=0.01, dc_type=:lar)...) - p_rf = slayer_parameters(; _ref_kwargs(; dr_val=0.01, dc_type=:rfitzp)...) - p_tor = slayer_parameters(; _ref_kwargs(; dr_val=0.01, dc_type=:toroidal)...) + p_lar = slayer_parameters(; _ref_kwargs(; delta_crit_D_R=0.01, delta_crit_type=:lar)...) + p_rf = slayer_parameters(; _ref_kwargs(; delta_crit_D_R=0.01, delta_crit_type=:fitzpatrick)...) + p_tor = slayer_parameters(; _ref_kwargs(; delta_crit_D_R=0.01, delta_crit_type=:toroidal)...) @test isfinite(p_lar.dc_tmp) @test isfinite(p_rf.dc_tmp) @test isfinite(p_tor.dc_tmp) - # dr_val > 0 with the (-dr_val) prefactor ⇒ negative dc_tmp for - # :lar, :rfitzp, :toroidal branches. + # delta_crit_D_R > 0 with the (-delta_crit_D_R) prefactor ⇒ negative dc_tmp for + # :lar, :fitzpatrick, :toroidal branches. @test p_lar.dc_tmp < 0 @test p_rf.dc_tmp < 0 @test p_tor.dc_tmp < 0 - # Sign flips with sign of dr_val + # Sign flips with sign of delta_crit_D_R p_lar_neg = slayer_parameters(; - _ref_kwargs(; dr_val=-0.01, dc_type=:lar)...) + _ref_kwargs(; delta_crit_D_R=-0.01, delta_crit_type=:lar)...) @test sign(p_lar_neg.dc_tmp) == -sign(p_lar.dc_tmp) - # Reject unknown dc_type + # Reject unknown delta_crit_type @test_throws ArgumentError slayer_parameters(; - _ref_kwargs(; dr_val=0.01, dc_type=:bogus)...) + _ref_kwargs(; delta_crit_D_R=0.01, delta_crit_type=:bogus)...) end @testset "Test 1c: SLAYERParameters direct kwarg construction" begin @@ -134,8 +134,8 @@ ) @test p.tau == 1.0 @test p.dc_tmp == 0.0 - @test p.dc_type === :none - @test p.dr_val == 0.0 + @test p.delta_crit_type === :none + @test p.delta_crit_D_R == 0.0 @test p.ising == 0 end diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index f74d774c9..9a1ba026d 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -37,7 +37,7 @@ @test_throws ArgumentError Runner.validate( SLAYERControl(; coupling_mode=:bogus)) @test_throws ArgumentError Runner.validate( - SLAYERControl(; dc_type=:bogus)) + SLAYERControl(; delta_crit_type=:bogus)) @test_throws ArgumentError Runner.validate( SLAYERControl(; msing_max=0)) @test_throws ArgumentError Runner.validate( @@ -50,11 +50,11 @@ "inner_model" => "slayer_fitzpatrick", "scan_mode" => "brute_force", "coupling_mode" => "coupled", - "dc_type" => "rfitzp", + "delta_crit_type" => "fitzpatrick", "msing_max" => 2, "bt" => 1.8, "mu_i" => 2.0, - "dr_val" => 0.01, + "delta_crit_D_R" => 0.01, "scan_grid" => Dict{String,Any}( "Q_re_range" => [-5.0, 5.0], "Q_im_range" => [-1.0, 3.0], @@ -73,10 +73,10 @@ @test c.inner_model === :slayer_fitzpatrick @test c.scan_mode === :brute_force @test c.coupling_mode === :coupled - @test c.dc_type === :rfitzp + @test c.delta_crit_type === :fitzpatrick @test c.msing_max == 2 @test c.bt === 1.8 - @test c.dr_val == 0.01 + @test c.delta_crit_D_R == 0.01 @test c.Q_re_range == (-5.0, 5.0) @test c.Q_im_range == (-1.0, 3.0) @test c.nre == 50 diff --git a/test/runtests_toml_backcompat.jl b/test/runtests_toml_backcompat.jl new file mode 100644 index 000000000..4ed15e3c1 --- /dev/null +++ b/test/runtests_toml_backcompat.jl @@ -0,0 +1,132 @@ +# Back-compat tests for deprecated TOML spellings: decks written with the old key/value +# names must load with a deprecation warning and produce control structs identical to the +# new spellings. +using Test +using TOML +using Logging +using GeneralizedPerturbedEquilibrium +using GeneralizedPerturbedEquilibrium.ForceFreeStates: ForceFreeStatesControl +using GeneralizedPerturbedEquilibrium.Equilibrium: EquilibriumConfig +using GeneralizedPerturbedEquilibrium.Runner: slayer_control_from_toml + +const GPE = GeneralizedPerturbedEquilibrium + +# Field-by-field struct equality (generic == is identity for mutable structs) +fields_equal(a::T, b::T) where {T} = all(getfield(a, f) == getfield(b, f) for f in fieldnames(T)) + +# Run f while discarding its log output (deprecation warnings are asserted separately) +quietly(f) = with_logger(f, NullLogger()) + +@testset "TOML deprecated-spelling back-compat" begin + @testset "_rename_keys! warns, remaps, and lets an explicit new key win" begin + t = Dict{String,Any}("newq0" => 1.5, "mpsi" => 64) + @test_logs (:warn, r"`newq0` in \[Equilibrium\] was renamed to `q0_override`") GPE._rename_keys!(t, GPE._RENAMED_EQUIL_KEYS, "Equilibrium") + @test !haskey(t, "newq0") + @test t["q0_override"] == 1.5 + @test t["mpsi"] == 64 + + t = Dict{String,Any}("newq0" => 1.5, "q0_override" => 2.5) + @test_logs (:warn, r"renamed to `q0_override`") GPE._rename_keys!(t, GPE._RENAMED_EQUIL_KEYS, "Equilibrium") + @test t["q0_override"] == 2.5 + end + + @testset "_rename_value! warns and remaps deprecated enum values" begin + t = Dict{String,Any}("f0type" => "jkp") + @test_logs (:warn, r"`f0type = \"jkp\"` in \[KineticForces\] is deprecated") GPE._rename_value!(t, "f0type", "jkp", "park", "KineticForces") + @test t["f0type"] == "park" + # Non-matching values pass through silently + t = Dict{String,Any}("f0type" => "park") + @test_logs GPE._rename_value!(t, "f0type", "jkp", "park", "KineticForces") + @test t["f0type"] == "park" + end + + @testset "integrator remap mirrors the old use_parallel/use_riccati dispatch" begin + # Old dispatch order: use_parallel (default true) wins, then use_riccati, else serial. + for (tbl, expect) in [ + Dict{String,Any}("use_parallel" => true) => "stride", + Dict{String,Any}("use_parallel" => false) => "serial", + Dict{String,Any}("use_parallel" => false, "use_riccati" => true) => "riccati", + Dict{String,Any}("use_riccati" => true) => "stride", + Dict{String,Any}("use_riccati" => false) => "stride" + ] + @test_logs (:warn, r"replaced by the `integrator` enum") GPE._remap_integrator_keys!(tbl) + @test tbl["integrator"] == expect + @test !haskey(tbl, "use_parallel") && !haskey(tbl, "use_riccati") + end + # An explicit integrator wins over the old flags + t = Dict{String,Any}("use_parallel" => false, "integrator" => "riccati") + @test_logs (:warn, r"ignored because `integrator` is also set") GPE._remap_integrator_keys!(t) + @test t["integrator"] == "riccati" + # No old keys: no warning, table untouched + t = Dict{String,Any}("integrator" => "stride") + @test_logs GPE._remap_integrator_keys!(t) + @test t["integrator"] == "stride" + end + + @testset "old FFS keys build an identical ForceFreeStatesControl" begin + old = Dict{String,Any}("use_parallel" => false, "parallel_threads" => 3, + "psiedge" => 0.97, "nstep" => 100, "diagnose_ca" => true, + "nn_low" => 1, "nn_high" => 1) + new = Dict{String,Any}("integrator" => "serial", "integrator_threads" => 3, + "dW_edge_scan_start" => 0.97, "nn_low" => 1, "nn_high" => 1) + quietly() do + GPE._rename_keys!(old, GPE._RENAMED_FFS_KEYS, "ForceFreeStates") + GPE._remap_integrator_keys!(old) + GPE._drop_deprecated_keys!(old, GPE._DEPRECATED_FFS_KEYS, "ForceFreeStates") + end + ctrl_old = ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in old)...) + ctrl_new = ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in new)...) + @test fields_equal(ctrl_old, ctrl_new) + end + + @testset "old Equilibrium keys/values build an identical EquilibriumConfig" begin + old = Dict{String,Any}("eq_type" => "efit", "eq_filename" => "g0.eqdsk", + "newq0" => 2, "use_galgrid" => false, "grid_type" => "ldp") + new = Dict{String,Any}("eq_type" => "efit", "eq_filename" => "g0.eqdsk", + "q0_override" => 2.0, "use_galerkin_grid" => false, "grid_type" => "rational_packed") + cfg_old = quietly() do + GPE._rename_keys!(old, GPE._RENAMED_EQUIL_KEYS, "Equilibrium") + EquilibriumConfig(old, ".") + end + cfg_new = quietly() do + EquilibriumConfig(new, ".") + end + @test cfg_old.q0_override == 2.0 + @test cfg_old.grid_type == "rational_packed" + @test fields_equal(cfg_old, cfg_new) + end + + @testset "old SLAYER keys/values build an identical SLAYERControl" begin + old = Dict{String,Any}("enabled" => true, "dc_type" => "rfitzp", + "dr_val" => 0.01, "dgeo_val" => 0.2) + new = Dict{String,Any}("enabled" => true, "delta_crit_type" => "fitzpatrick", + "delta_crit_D_R" => 0.01, "delta_crit_geo_factor" => 0.2) + ctrl_old = quietly() do + slayer_control_from_toml(old) + end + ctrl_new = slayer_control_from_toml(new) + @test ctrl_old.delta_crit_type === :fitzpatrick + @test fields_equal(ctrl_old, ctrl_new) + # The rename warnings actually fire + @test_logs (:warn, r"`dc_type` in \[SLAYER\] was renamed") match_mode = :any slayer_control_from_toml(Dict{String,Any}("dc_type" => "lar")) + end + + @testset "build_inputs_from_toml applies the Equilibrium renames on a real deck" begin + mktempdir() do dir + write(joinpath(dir, "gpec.toml"), + """ + [Equilibrium] + eq_type = "efit" + eq_filename = "g_unused.eqdsk" + newq0 = 0 + use_galgrid = true + """) + inputs, eq_config, _ = quietly() do + GPE.build_inputs_from_toml(dir) + end + @test !haskey(inputs["Equilibrium"], "newq0") + @test eq_config.q0_override == 0.0 + @test eq_config.use_galerkin_grid === true + end + end +end diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 630c60436..578dd6e1e 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -465,7 +465,7 @@ @testset "extract_plasma_surface_at_psi" begin # Self-contained analytic Solovev equilibrium (same recipe as runtests_equil.jl). - eq_config = Equilibrium.EquilibriumConfig(; eq_type="sol", eq_filename="unused", jac_type="pest", grid_type="ldp", psilow=1e-4, psihigh=0.99999, mpsi=64, mtheta=128) + eq_config = Equilibrium.EquilibriumConfig(; eq_type="sol", eq_filename="unused", jac_type="pest", grid_type="rational_packed", psilow=1e-4, psihigh=0.99999, mpsi=64, mtheta=128) sol_config = Equilibrium.SolovevConfig(64, 64, 64, 1.6, 0.33, 1.0, 1.9, 1.0, 1.0, 1.0) pe = Equilibrium.equilibrium_solver(Equilibrium.sol_run(eq_config, sol_config)) diff --git a/test/test_data/regression_solovev_ideal_example/gpec.toml b/test/test_data/regression_solovev_ideal_example/gpec.toml index a16ab6898..28e5fb9fe 100644 --- a/test/test_data/regression_solovev_ideal_example/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example/gpec.toml @@ -5,12 +5,12 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 16 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -28,7 +28,7 @@ equal_arc_wall = true # Equal arc length distribution of nodes local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) @@ -46,8 +46,8 @@ 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) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (unlocks SingularSurfaces/delta_prime_matrix) +integrator_threads = 2 # Stride 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 dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml index 0f0bc5c47..3e939e62a 100644 --- a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml @@ -5,12 +5,12 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 16 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -28,7 +28,7 @@ equal_arc_wall = true # Equal arc length distribution of nodes local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) @@ -46,8 +46,8 @@ 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 (multi-n Δ' matrix has open issues — sing_lim! warns and skips — but ξ and energies are valid) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (multi-n Δ' matrix has open issues — sing_lim! warns and skips — but ξ and energies are valid) +integrator_threads = 2 # Stride 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 multi-n — dmlim truncation is ambiguous when n varies (sing_lim! skips anyway) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml index 5b87267a9..f8df5d802 100644 --- a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml @@ -5,12 +5,12 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 16 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -29,7 +29,7 @@ equal_arc_wall = true # Equal arc length distribution of nodes local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) diff --git a/test/test_data/regression_solovev_kinetic_example/gpec.toml b/test/test_data/regression_solovev_kinetic_example/gpec.toml index 1c020a086..9f7235e86 100644 --- a/test/test_data/regression_solovev_kinetic_example/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_example/gpec.toml @@ -5,12 +5,12 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 16 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -28,7 +28,7 @@ equal_arc_wall = true # Equal arc length distribution of nodes local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) @@ -46,8 +46,8 @@ 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) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (unlocks SingularSurfaces/delta_prime_matrix) +integrator_threads = 2 # Stride 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 dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml index b11caab8b..c80cc8aa9 100644 --- a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml @@ -5,12 +5,12 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 16 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -28,7 +28,7 @@ equal_arc_wall = true # Equal arc length distribution of nodes local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) @@ -46,8 +46,8 @@ 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 (multi-n Δ' matrix has open issues — sing_lim! warns and skips — but ξ and energies are valid) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) +integrator = "stride" # Integration algorithm: "stride" FM-propagator BVP (multi-n Δ' matrix has open issues — sing_lim! warns and skips — but ξ and energies are valid) +integrator_threads = 2 # Stride 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 multi-n — dmlim truncation is ambiguous when n varies (sing_lim! skips anyway) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml index ca6eeb003..825e561b0 100644 --- a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml @@ -6,12 +6,12 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "ldp" # Radial grid packing type +grid_type = "rational_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized flux coordinate psihigh = 0.9995 # Upper limit of normalized flux coordinate mpsi = 16 # Number of radial grid points mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +q0_override = 0.0 # Override for on-axis safety factor (0 = use input value) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) @@ -29,7 +29,7 @@ equal_arc_wall = true # Equal arc length distribution of nodes local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -psiedge = 0.99 # Edge dW scan band: dW(ψ) computed for ψ ∈ [psiedge, psilim], integration truncated at peak +dW_edge_scan_start = 0.99 # Edge dW scan band: dW(ψ) computed for ψ ∈ [dW_edge_scan_start, psilim], integration truncated at peak qlow = 1.02 # Integration initiated at q determined by min(q0, qlow)... qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh)... sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) From cab894a0052ca5586ee777ba11531d8f3dfeebe4 Mon Sep 17 00:00:00 2001 From: priyanshlunia <40486607+priyanshlunia@users.noreply.github.com> Date: Fri, 14 Aug 2026 16:47:10 -0400 Subject: [PATCH 4/4] FORCINGTERMS - REFACTOR - Rename coil_set xnom/ynom/znom to rotation_center_x/_y/_z (issue #287) Completes the issue #287 rename sweep: the [[ForcingTerms.coil_set]] explicit rotation-center keys move to descriptive spellings, with the same warn-and-remap back-compat as the other renames (the coil parser previously dropped unknown keys silently, so the remap is required for old decks to keep working). Old keys added to the pre-commit toml-no-deprecated-keys pattern and covered in runtests_toml_backcompat.jl. | Section | Old | New | |----------------------------|------|-------------------| | [[ForcingTerms.coil_set]] | xnom | rotation_center_x | | [[ForcingTerms.coil_set]] | ynom | rotation_center_y | | [[ForcingTerms.coil_set]] | znom | rotation_center_z | Co-Authored-By: Claude Fable 5 --- .pre-commit-config.yaml | 2 +- src/ForcingTerms/CoilGeometry.jl | 35 +++++++++++++++++++------------- test/runtests_toml_backcompat.jl | 17 ++++++++++++++-- 3 files changed, 37 insertions(+), 17 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 506dd06ff..94c59f489 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -69,5 +69,5 @@ repos: - id: toml-no-deprecated-keys name: 'TOML conventions: no deprecated config keys' language: pygrep - entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|power_bp|power_b|power_r|power_rc|nstep|diagnose_ca|use_parallel|use_riccati|parallel_threads|psiedge|newq0|use_galgrid|reg_spot|dc_type|dr_val|dgeo_val)\s*=' + entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|power_bp|power_b|power_r|power_rc|nstep|diagnose_ca|use_parallel|use_riccati|parallel_threads|psiedge|newq0|use_galgrid|reg_spot|dc_type|dr_val|dgeo_val|xnom|ynom|znom)\s*=' files: ^(examples/.*\.toml|test/test_data/.*\.toml)$ diff --git a/src/ForcingTerms/CoilGeometry.jl b/src/ForcingTerms/CoilGeometry.jl index b2d19ab21..cd5824069 100644 --- a/src/ForcingTerms/CoilGeometry.jl +++ b/src/ForcingTerms/CoilGeometry.jl @@ -69,7 +69,7 @@ Shifts/tilts apply to every source (the analytic geometry is built first, then t - `currents`: current [A] per conductor `[ncoil]`; shorter arrays pad with zeros - `shiftx`, `shifty`, `shiftz`: per-conductor translation [m] `[ncoil]` - `tiltx`, `tilty`, `tiltz`: per-conductor tilt in degrees (or meters if `tilt_in_meters`) - - `xnom`, `ynom`, `znom`: explicit rotation center [m]; defaults to arc-length-weighted center of mass + - `rotation_center_x`, `rotation_center_y`, `rotation_center_z`: explicit rotation center [m]; defaults to arc-length-weighted center of mass - `n_tilt`: toroidal mode number for tilt/shift modulation; -1 means inherit run's n - `tilt_in_meters`: interpret tilt as displacement [m] instead of angle [degrees] @@ -106,9 +106,9 @@ Base.@kwdef struct CoilSetConfig tiltx::Vector{Float64} = Float64[] tilty::Vector{Float64} = Float64[] tiltz::Vector{Float64} = Float64[] - xnom::Vector{Float64} = Float64[] - ynom::Vector{Float64} = Float64[] - znom::Vector{Float64} = Float64[] + rotation_center_x::Vector{Float64} = Float64[] + rotation_center_y::Vector{Float64} = Float64[] + rotation_center_z::Vector{Float64} = Float64[] n_tilt::Int = -1 tilt_in_meters::Bool = false @@ -180,6 +180,13 @@ function CoilConfig(ft_ctrl::ForcingTermsControl) end function _parse_coil_set_config(d::Dict{String,Any}) + # Deprecated key spellings, accepted with a warning until removal after v2.0.0. + for (old, new) in ("xnom" => "rotation_center_x", "ynom" => "rotation_center_y", "znom" => "rotation_center_z") + haskey(d, old) || continue + @warn "`$old` in [[ForcingTerms.coil_set]] was renamed to `$new`; the old key is deprecated and will be removed after v2.0.0." + haskey(d, new) || (d[new] = d[old]) + delete!(d, old) + end fvec(key) = Float64.(get(d, key, Float64[])) # rz_corners arrives as a Vector of [R, Z] pairs (TOML array of arrays) corners = [Float64.(c) for c in get(d, "rz_corners", Vector{Float64}[])] @@ -195,9 +202,9 @@ function _parse_coil_set_config(d::Dict{String,Any}) tiltx=fvec("tiltx"), tilty=fvec("tilty"), tiltz=fvec("tiltz"), - xnom=fvec("xnom"), - ynom=fvec("ynom"), - znom=fvec("znom"), + rotation_center_x=fvec("rotation_center_x"), + rotation_center_y=fvec("rotation_center_y"), + rotation_center_z=fvec("rotation_center_z"), n_tilt=get(d, "n_tilt", -1), tilt_in_meters=get(d, "tilt_in_meters", false), radius=Float64(get(d, "radius", 0.0)), @@ -777,7 +784,7 @@ Apply per-conductor shifts and tilts to a coil set, returning a modified copy. Replicates the Fortran `coil_read` shift/tilt logic (coil.F lines 240–340): - - Tilts are rotations around the arc-length-weighted center of mass (unless `xnom/ynom/znom` specified) + - Tilts are rotations around the arc-length-weighted center of mass (unless `rotation_center_x/rotation_center_y/rotation_center_z` specified) - `n_tilt` controls the toroidal periodicity of tilt/shift modulation - n_tilt = 0: rigid shift only (no tilts applied) - n_tilt ≥ 1: n-fold modulated perturbations @@ -797,9 +804,9 @@ function apply_transforms(cs::CoilSet, cfg::CoilSetConfig; n_tilt::Int=1) tilty_cfg = _pad(cfg.tilty, ncoil) tiltz_cfg = _pad(cfg.tiltz, ncoil) - xnom_cfg = _pad(isempty(cfg.xnom) ? fill(_NOM_UNSET_SENTINEL, ncoil) : cfg.xnom, ncoil) - ynom_cfg = _pad(isempty(cfg.ynom) ? fill(_NOM_UNSET_SENTINEL, ncoil) : cfg.ynom, ncoil) - znom_cfg = _pad(isempty(cfg.znom) ? fill(_NOM_UNSET_SENTINEL, ncoil) : cfg.znom, ncoil) + rotation_center_x_cfg = _pad(isempty(cfg.rotation_center_x) ? fill(_NOM_UNSET_SENTINEL, ncoil) : cfg.rotation_center_x, ncoil) + rotation_center_y_cfg = _pad(isempty(cfg.rotation_center_y) ? fill(_NOM_UNSET_SENTINEL, ncoil) : cfg.rotation_center_y, ncoil) + rotation_center_z_cfg = _pad(isempty(cfg.rotation_center_z) ? fill(_NOM_UNSET_SENTINEL, ncoil) : cfg.rotation_center_z, ncoil) # Check if n_tilt = 0 suppresses tilts (Fortran: "no n=0 component") apply_tilt = n_tilt != 0 @@ -827,9 +834,9 @@ function apply_transforms(cs::CoilSet, cfg::CoilSetConfig; n_tilt::Int=1) view(cs.x, j, k, :), view(cs.y, j, k, :), view(cs.z, j, k, :) ) # Use user-specified center if provided (|nom| < _NOM_THRESHOLD, matching Fortran) - x0 = abs(xnom_cfg[j]) < _NOM_THRESHOLD ? xnom_cfg[j] : cx - y0 = abs(ynom_cfg[j]) < _NOM_THRESHOLD ? ynom_cfg[j] : cy - z0 = abs(znom_cfg[j]) < _NOM_THRESHOLD ? znom_cfg[j] : cz + x0 = abs(rotation_center_x_cfg[j]) < _NOM_THRESHOLD ? rotation_center_x_cfg[j] : cx + y0 = abs(rotation_center_y_cfg[j]) < _NOM_THRESHOLD ? rotation_center_y_cfg[j] : cy + z0 = abs(rotation_center_z_cfg[j]) < _NOM_THRESHOLD ? rotation_center_z_cfg[j] : cz (x0, y0, z0) end diff --git a/test/runtests_toml_backcompat.jl b/test/runtests_toml_backcompat.jl index 4ed15e3c1..44699c35b 100644 --- a/test/runtests_toml_backcompat.jl +++ b/test/runtests_toml_backcompat.jl @@ -11,8 +11,8 @@ using GeneralizedPerturbedEquilibrium.Runner: slayer_control_from_toml const GPE = GeneralizedPerturbedEquilibrium -# Field-by-field struct equality (generic == is identity for mutable structs) -fields_equal(a::T, b::T) where {T} = all(getfield(a, f) == getfield(b, f) for f in fieldnames(T)) +# Field-by-field struct equality (generic == is identity for mutable structs; isequal so NaN sentinels compare equal) +fields_equal(a::T, b::T) where {T} = all(isequal(getfield(a, f), getfield(b, f)) for f in fieldnames(T)) # Run f while discarding its log output (deprecation warnings are asserted separately) quietly(f) = with_logger(f, NullLogger()) @@ -111,6 +111,19 @@ quietly(f) = with_logger(f, NullLogger()) @test_logs (:warn, r"`dc_type` in \[SLAYER\] was renamed") match_mode = :any slayer_control_from_toml(Dict{String,Any}("dc_type" => "lar")) end + @testset "old coil_set keys build an identical CoilSetConfig" begin + old = Dict{String,Any}("name" => "c79", "xnom" => [1.0], "ynom" => [2.0], "znom" => [3.0]) + new = Dict{String,Any}("name" => "c79", "rotation_center_x" => [1.0], + "rotation_center_y" => [2.0], "rotation_center_z" => [3.0]) + cfg_old = quietly() do + GPE.ForcingTerms._parse_coil_set_config(old) + end + cfg_new = GPE.ForcingTerms._parse_coil_set_config(new) + @test cfg_old.rotation_center_x == [1.0] + @test fields_equal(cfg_old, cfg_new) + @test_logs (:warn, r"`xnom` in \[\[ForcingTerms.coil_set\]\] was renamed") match_mode = :any GPE.ForcingTerms._parse_coil_set_config(Dict{String,Any}("xnom" => [1.0])) + end + @testset "build_inputs_from_toml applies the Equilibrium renames on a real deck" begin mktempdir() do dir write(joinpath(dir, "gpec.toml"),