diff --git a/CLAUDE.md b/CLAUDE.md index cbd8e2cad..206b6e39b 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -123,6 +123,14 @@ Roster, the recommended review pipeline (physics fidelity → readability → pe ### Plotting Spectrum-plot (`seriestype=:steppre` / `step_series`) and figure-saving conventions are in **[`docs/development/plotting.md`](docs/development/plotting.md)** — read before writing any plotting code. +**Adding a plot.** Every plot in `src/Analysis/` takes a result in memory or a `gpec.h5` path and is registered in `src/Analysis/Registry.jl`, which `available_plots` and the REPL summaries read: +- Read the data through one `_x_data(src)` method per source (the result's fields, and the same datasets from the file), and draw from the returned NamedTuple only, so both sources go through one plotting body. +- Start the public `plot_x(src; ...)` with `_require(plot_x, src)`, which throws the registry's reason when the source lacks the data; never draw a "No data" placeholder or return `nothing`. Composite summaries skip the panels they cannot draw (`_can`) and throw only when none can. +- Call `register!(plot_x, sources, needs, description)` after the definition. `needs` names capabilities; a new one gets a `capabilities` predicate for each result type that can hold it, a dataset in `_H5_CAPABILITIES`, and a `REASONS` sentence saying what to change in the run. +- `test/runtests_analysis_plots.jl` checks every registered plot from both sources automatically; add an extractor to its list when the plot reads new data. + +**REPL summaries.** `src/Display.jl` gives every result type a one-line `show` and a short `text/plain` summary. A new result type gets both, plus an entry in `_DISPLAY_FIELDS`; a new field on a summarized type is classified there as shown or omitted, which `test/display_check.jl` enforces. + ### Subagent Consultations When delegating to specialized agents (julia-performance-optimizer, fast-interpolations-optimizer, fortran-physics-reviewer, clean-code-reviewer, etc.), every prompt **must include an explicit budget** because runaway agents silently consume the user's daily token quota. A single consultation that explores instead of editing has been measured at 167 tool calls / 55 minutes / no return — that is a session-killer. Defaults: diff --git a/docs/src/analysis.md b/docs/src/analysis.md index 6034f2caf..9ca4ceb88 100644 --- a/docs/src/analysis.md +++ b/docs/src/analysis.md @@ -1,6 +1,6 @@ # Analysis -`Analysis` holds the plotting and post-processing functions for GPEC results, grouped by the stage they read; each works on the in-memory result or on `gpec.h5`. +`Analysis` holds the plotting and post-processing functions for GPEC results, grouped by the stage they read. Each registered plot takes the in-memory result or the path of the `gpec.h5` it was saved to; `available_plots` lists them with what each source can draw (see [Plotting a result](api.md#Plotting-a-result)). ## Submodules @@ -18,6 +18,13 @@ Figure conventions for spectrum plots are in `docs/development/plotting.md` of t GeneralizedPerturbedEquilibrium.Analysis ``` +## Plot registry + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Analysis] +Pages = ["Registry.jl"] +``` + ## ForceFreeStates ```@autodocs diff --git a/docs/src/api.md b/docs/src/api.md index 06a117a30..d12110902 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -53,10 +53,31 @@ ForceFreeStatesResult Riccati, ideal closure Δ′_ii (ψ_N) 2/1: 0.1781, 3/1: -215.7 Solution not computed (Riccati stores no ξ profiles) Wall ConformalWall, standoff = 0.2415 a + Plots 11 of 13 available (see available_plots) ``` The fields themselves are documented in each type's docstring. +## Plotting a result + +Every plot in [`Analysis`](analysis.md) takes the result in memory or the path of the `gpec.h5` it was saved to. `available_plots` lists the plots for a result, a path, or the named tuple `main` returns, and for each one the run cannot draw, what to change: + +```julia-repl +julia> available_plots(ffs) +Plots for ForceFreeStatesResult (riccati): 11 of 13 available + ✓ ForceFreeStates.plot_ffs_summary the stability panels the run supplies + ✗ ForceFreeStates.plot_mode_displacement needs stored ξ profiles: use the Forward integrator or a matched Galerkin solve (Riccati and unmatched Galerkin store none) + ✓ ForceFreeStates.plot_fixed_boundary_stability_criterion critical determinant against ψ_N + ⋮ +``` + +A plot the source cannot draw throws an `ArgumentError` with the same reason, and a summary leaves out the panels it cannot draw: + +```julia +Analysis.ForceFreeStates.plot_ffs_summary(ffs) # the result in memory +Analysis.ForceFreeStates.plot_ffs_summary("gpec.h5") # or the file to which it was saved +``` + ## Saving and rerunning `save_to_hdf5(target, result)` saves any stage result; `target` is a path, which opens or creates the file, or an open HDF5 file or group, which stays owned by the caller. diff --git a/examples/Solovev_ideal_example/analyze_example.jl b/examples/Solovev_ideal_example/analyze_example.jl index 62178b316..91400a204 100644 --- a/examples/Solovev_ideal_example/analyze_example.jl +++ b/examples/Solovev_ideal_example/analyze_example.jl @@ -4,15 +4,14 @@ using GeneralizedPerturbedEquilibrium, Plots using GeneralizedPerturbedEquilibrium: Analysis isinteractive() ? plotlyjs() : gr() -h5path = joinpath(@__DIR__, "gpec.h5") +# Run the deck and plot its results in memory; every plot also takes the path of the gpec.h5 the run wrote. +ret = GeneralizedPerturbedEquilibrium.main([@__DIR__]) +display(available_plots(ret)) -# Summary plots -p_eq = Analysis.Equilibrium.plot_equilibrium_summary(h5path) -p_ffs = Analysis.ForceFreeStates.plot_ffs_summary(h5path) +p_eq = Analysis.Equilibrium.plot_equilibrium_summary(ret.equil) +p_ffs = Analysis.ForceFreeStates.plot_ffs_summary(ret.ffs) display(p_eq); Plots.savefig(p_eq, joinpath(@__DIR__, "equilibrium_summary.png")); display(p_ffs); Plots.savefig(p_ffs, joinpath(@__DIR__, "ffs_summary.png")); -display(p_pe); -Plots.savefig(p_pe, joinpath(@__DIR__, "pe_summary.png")); diff --git a/src/Analysis/Analysis.jl b/src/Analysis/Analysis.jl index 59309953d..cf3e8abdc 100644 --- a/src/Analysis/Analysis.jl +++ b/src/Analysis/Analysis.jl @@ -3,6 +3,10 @@ Post-processing and visualization utilities for GPEC simulation outputs. +Every registered plot takes a result in memory or the path of the `gpec.h5` it was saved to, and throws an +`ArgumentError` naming what to change in the run when that source lacks the data it draws. `available_plots(src)` +lists the plots for a result, a path, or the named tuple `main` returns. + ## Submodules - `ForceFreeStates`: Plotting functions for ForceFreeStates (DCON-style ideal MHD stability) results @@ -14,6 +18,7 @@ Post-processing and visualization utilities for GPEC simulation outputs. """ module Analysis +include("Registry.jl") include("ForceFreeStates.jl") include("Equilibrium.jl") include("CoilForcing.jl") @@ -21,4 +26,6 @@ include("PerturbedEquilibrium.jl") include("PerturbedEquilibriumModes.jl") include("ErrorFields.jl") +export available_plots + end # module Analysis diff --git a/src/Analysis/CoilForcing.jl b/src/Analysis/CoilForcing.jl index 660d0505b..371c5548f 100644 --- a/src/Analysis/CoilForcing.jl +++ b/src/Analysis/CoilForcing.jl @@ -17,6 +17,8 @@ module CoilForcing using Plots import ...Utilities.FourierTransforms: FourierTransform, inverse +import ...PerturbedEquilibrium: PerturbedEquilibriumState +using ..Analysis: register!, _require """ plot_coil_geometry_3d(coil_sets; save_path=nothing, kwargs...) @@ -453,6 +455,33 @@ function plot_surface_3d(equil, modes, n::Int; coil_sets=nothing, mtheta::Int=18 psi=psi, colorbar_title=colorbar_title, save_path=save_path, kwargs...) end +""" + plot_coil_geometry_3d(pe::PerturbedEquilibriumState; kwargs...) + plot_coil_geometry_rz(pe::PerturbedEquilibriumState; kwargs...) + plot_mode_spectrum(pe::PerturbedEquilibriumState; kwargs...) + +The coil plots of the coil sets a perturbed equilibrium was forced by: the (R, Z) cross-sections over its plasma +boundary, and the forcing spectrum across its poloidal mode range. +""" +function plot_coil_geometry_3d(pe::PerturbedEquilibriumState; kwargs...) + _require(plot_coil_geometry_3d, pe) + return plot_coil_geometry_3d(pe.coil_sets; kwargs...) +end + +function plot_coil_geometry_rz(pe::PerturbedEquilibriumState; kwargs...) + _require(plot_coil_geometry_rz, pe) + return plot_coil_geometry_rz(pe.coil_sets; equil=pe.response.ffs.equil, kwargs...) +end + +function plot_mode_spectrum(pe::PerturbedEquilibriumState; kwargs...) + _require(plot_mode_spectrum, pe) + return plot_mode_spectrum(pe.forcing_modes; mlow=pe.response.ffs.mlow, mhigh=pe.response.ffs.mhigh, kwargs...) +end + +register!(plot_coil_geometry_3d, (PerturbedEquilibriumState,), (:coil_sets,), "3D coil conductors") +register!(plot_coil_geometry_rz, (PerturbedEquilibriumState,), (:coil_sets,), "coil cross-sections over the plasma boundary") +register!(plot_mode_spectrum, (PerturbedEquilibriumState,), (:forcing_modes,), "|Φ_x(m, n)| forcing spectrum") + export plot_coil_geometry_3d, plot_coil_geometry_rz export plot_bn_contour, plot_mode_spectrum export plot_surface_3d, control_surface_scalar diff --git a/src/Analysis/Equilibrium.jl b/src/Analysis/Equilibrium.jl index f00a406cd..699158851 100644 --- a/src/Analysis/Equilibrium.jl +++ b/src/Analysis/Equilibrium.jl @@ -1,7 +1,8 @@ """ Equilibrium -Post-processing and visualization functions for GPEC equilibrium objects and HDF5 outputs. +Post-processing and visualization functions for GPEC equilibria. Each plot takes a `PlasmaEquilibrium`, a +`ForceFreeStatesResult` (which adds its rational surfaces as markers), or the path of a `gpec.h5`. """ module Equilibrium @@ -9,15 +10,84 @@ using HDF5 using LaTeXStrings using Plots +using ..Analysis: register!, _require, _can, capabilities +import ...Equilibrium: PlasmaEquilibrium +import ...ForceFreeStates: ForceFreeStatesResult + +const Source = Union{AbstractString,PlasmaEquilibrium,ForceFreeStatesResult} + +# Rational surfaces (ψ_N, q) to mark: none for a bare equilibrium. +_surfaces(::PlasmaEquilibrium) = (; sing_psi=Float64[], sing_q=Float64[]) +_surfaces(r::ForceFreeStatesResult) = (; sing_psi=[s.psifac for s in r.surfaces], sing_q=[s.q for s in r.surfaces]) +function _surfaces(fid::HDF5.File) + haskey(fid, "SingularSurfaces/rational_psi") || return _surfaces(nothing) + return (; sing_psi=read(fid["SingularSurfaces/rational_psi"]), sing_q=read(fid["SingularSurfaces/rational_q"])) +end +_surfaces(::Nothing) = (; sing_psi=Float64[], sing_q=Float64[]) + +_equil(eq::PlasmaEquilibrium) = eq +_equil(r::ForceFreeStatesResult) = r.equil + +function _profile_data(src::Union{PlasmaEquilibrium,ForceFreeStatesResult}) + eq = _equil(src) + prof = eq.profiles + return (; psi=prof.xs, q=prof.q_spline.y, mu0p=prof.P_spline.y, twopif=prof.F_spline.y, q0=eq.params.q0, q95=eq.params.q95, _surfaces(src)...) +end + +function _profile_data(h5path::AbstractString) + h5open(h5path, "r") do fid + (; psi=read(fid["Equilibrium/Profiles/psi"]), q=read(fid["Equilibrium/Profiles/q"]), mu0p=read(fid["Equilibrium/Profiles/mu0p"]), + twopif=read(fid["Equilibrium/Profiles/2piF"]), q0=read(fid["Equilibrium/q_axis"]), q95=read(fid["Equilibrium/q_95"]), _surfaces(fid)...) + end +end + +function _geometry_data(src::Union{PlasmaEquilibrium,ForceFreeStatesResult}) + eq = _equil(src) + return (; rcoords=eq.rzphi_rsquared.nodal_derivs.partials[1, :, :], offset=eq.rzphi_offset.nodal_derivs.partials[1, :, :], + psi=eq.rzphi_xs, theta=eq.rzphi_ys, ro=eq.ro, zo=eq.zo, _surfaces(src)...) +end + +function _geometry_data(h5path::AbstractString) + h5open(h5path, "r") do fid + (; rcoords=read(fid["Equilibrium/Geometry/rcoords"]), offset=read(fid["Equilibrium/Geometry/offset"]), psi=read(fid["Equilibrium/Geometry/psi"]), + theta=read(fid["Equilibrium/Geometry/theta"]), ro=read(fid["Equilibrium/R_axis"]), zo=read(fid["Equilibrium/Z_axis"]), _surfaces(fid)...) + end +end + +# R, Z on the full nodal (ψ, θ) grid of the geometry data. +function _rz_grid(d) + R = Matrix{Float64}(undef, length(d.psi), length(d.theta)) + Z = similar(R) + for ipsi in eachindex(d.psi), itheta in eachindex(d.theta) + rfac = sqrt(max(0.0, d.rcoords[ipsi, itheta])) + eta = 2π * (d.theta[itheta] + d.offset[ipsi, itheta]) + R[ipsi, itheta] = d.ro + rfac * cos(eta) + Z[ipsi, itheta] = d.zo + rfac * sin(eta) + end + return R, Z +end + +function _global_data(src::Union{PlasmaEquilibrium,ForceFreeStatesResult}) + p = _equil(src).params + return (; q0=p.q0, q95=p.q95, betat=p.betat, betan=p.betan, kappa=p.kappa, li1=p.li1) +end + +function _global_data(h5path::AbstractString) + h5open(h5path, "r") do fid + (; q0=read(fid["Equilibrium/q_axis"]), q95=read(fid["Equilibrium/q_95"]), betat=read(fid["Equilibrium/beta_t"]), + betan=read(fid["Equilibrium/beta_N"]), kappa=read(fid["Equilibrium/kappa"]), li1=read(fid["Equilibrium/l_i_1"])) + end +end + """ - plot_qprofile(h5path; show_singular=true, save_path=nothing) + plot_qprofile(src; show_singular=true, save_path=nothing) Plot the safety factor q(ψ) profile, with optional vertical markers at each rational surface and horizontal reference lines at q0 and q95. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file + - `src`: A `PlasmaEquilibrium`, a `ForceFreeStatesResult`, or the path of a `gpec.h5` ### Keyword arguments @@ -28,12 +98,15 @@ and horizontal reference lines at q0 and q95. A `Plots.jl` plot object. """ -function plot_qprofile(h5path; show_singular=true, save_path=nothing) - xs, q, q0, q95 = h5open(h5path, "r") do fid - read(fid["Equilibrium/Profiles/psi"]), read(fid["Equilibrium/Profiles/q"]), - read(fid["Equilibrium/q_axis"]), read(fid["Equilibrium/q_95"]) - end +function plot_qprofile(src::Source; show_singular=true, save_path=nothing) + _require(plot_qprofile, src) + p = _plot_qprofile(_profile_data(src); show_singular) + isnothing(save_path) || savefig(p, save_path) + return p +end +function _plot_qprofile(d; show_singular=true) + xs, q, q0, q95 = d.psi, d.q, d.q0, d.q95 p = plot( xs, q; xlabel="Norm. Poloidal Flux", @@ -50,28 +123,23 @@ function plot_qprofile(h5path; show_singular=true, save_path=nothing) annotate!(p, 1.0, q95, text(" q95=$(round(q95, digits=2))", 7, :left, :gray)) if show_singular - msing, psi_sing, q_sing = h5open(h5path, "r") do fid - read(fid["SingularSurfaces/rational_count"]), read(fid["SingularSurfaces/rational_psi"]), read(fid["SingularSurfaces/rational_q"]) - end - for s in 1:msing - vline!(p, [psi_sing[s]]; linestyle=:dash, color=:red, label=nothing) - annotate!(p, psi_sing[s], q_sing[s], - text(" q=$(round(q_sing[s], digits=2))", 7, :left, :red)) + for (psi_s, q_s) in zip(d.sing_psi, d.sing_q) + vline!(p, [psi_s]; linestyle=:dash, color=:red, label=nothing) + annotate!(p, psi_s, q_s, + text(" q=$(round(q_s, digits=2))", 7, :left, :red)) end end - - isnothing(save_path) || savefig(p, save_path) return p end """ - plot_pressure_profile(h5path; save_path=nothing) + plot_pressure_profile(src; save_path=nothing) Plot the μ₀p(ψ) pressure profile. Vertical dashed lines mark rational surfaces if present. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file + - `src`: A `PlasmaEquilibrium`, a `ForceFreeStatesResult`, or the path of a `gpec.h5` ### Keyword arguments @@ -81,14 +149,16 @@ Plot the μ₀p(ψ) pressure profile. Vertical dashed lines mark rational surfac A `Plots.jl` plot object. """ -function plot_pressure_profile(h5path; save_path=nothing) - xs, mu0p, msing, psi_sing = h5open(h5path, "r") do fid - read(fid["Equilibrium/Profiles/psi"]), read(fid["Equilibrium/Profiles/mu0p"]), - read(fid["SingularSurfaces/rational_count"]), read(fid["SingularSurfaces/rational_psi"]) - end +function plot_pressure_profile(src::Source; save_path=nothing) + _require(plot_pressure_profile, src) + p = _plot_pressure_profile(_profile_data(src)) + isnothing(save_path) || savefig(p, save_path) + return p +end +function _plot_pressure_profile(d) p = plot( - xs, mu0p; + d.psi, d.mu0p; xlabel="Norm. Poloidal Flux", ylabel="μ₀p", title="", @@ -97,22 +167,20 @@ function plot_pressure_profile(h5path; save_path=nothing) left_margin=10Plots.mm, bottom_margin=5Plots.mm ) - for s in 1:msing - vline!(p, [psi_sing[s]]; linestyle=:dash, color=:red, label=nothing) + for psi_s in d.sing_psi + vline!(p, [psi_s]; linestyle=:dash, color=:red, label=nothing) end - - isnothing(save_path) || savefig(p, save_path) return p end """ - plot_f_profile(h5path; save_path=nothing) + plot_f_profile(src; save_path=nothing) Plot the toroidal field function 2πF(ψ) profile (F = RBφ/(2π)). ### Arguments - - `h5path`: Path to a GPEC HDF5 output file + - `src`: A `PlasmaEquilibrium`, a `ForceFreeStatesResult`, or the path of a `gpec.h5` ### Keyword arguments @@ -122,14 +190,16 @@ Plot the toroidal field function 2πF(ψ) profile (F = RBφ/(2π)). A `Plots.jl` plot object. """ -function plot_f_profile(h5path; save_path=nothing) - xs, twopif, msing, psi_sing = h5open(h5path, "r") do fid - read(fid["Equilibrium/Profiles/psi"]), read(fid["Equilibrium/Profiles/2piF"]), - read(fid["SingularSurfaces/rational_count"]), read(fid["SingularSurfaces/rational_psi"]) - end +function plot_f_profile(src::Source; save_path=nothing) + _require(plot_f_profile, src) + p = _plot_f_profile(_profile_data(src)) + isnothing(save_path) || savefig(p, save_path) + return p +end +function _plot_f_profile(d) p = plot( - xs, twopif; + d.psi, d.twopif; xlabel="Norm. Poloidal Flux", ylabel="2πF", title="", @@ -138,26 +208,24 @@ function plot_f_profile(h5path; save_path=nothing) left_margin=10Plots.mm, bottom_margin=5Plots.mm ) - for s in 1:msing - vline!(p, [psi_sing[s]]; linestyle=:dash, color=:red, label=nothing) + for psi_s in d.sing_psi + vline!(p, [psi_s]; linestyle=:dash, color=:red, label=nothing) end - - isnothing(save_path) || savefig(p, save_path) return p end """ - plot_flux_surfaces(h5path; n_psi=11, n_theta=18, save_path=nothing) + plot_flux_surfaces(src; n_psi=11, n_theta=18, save_path=nothing) Plot flux surface contours (constant ψ, blue) and field-line angle spokes (constant θ, red) -in physical (R, Z) space, reading nodal grid data directly from HDF5. +in physical (R, Z) space, from the nodal grid data. Psi contours are drawn at `n_psi` evenly spaced values between psilow and psihigh. Theta spokes are drawn at `n_theta` evenly spaced values. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file + - `src`: A `PlasmaEquilibrium`, a `ForceFreeStatesResult`, or the path of a `gpec.h5` ### Keyword arguments @@ -169,28 +237,19 @@ Theta spokes are drawn at `n_theta` evenly spaced values. A `Plots.jl` plot object. """ -function plot_flux_surfaces(h5path; n_psi=11, n_theta=18, save_path=nothing) - rcoords, offset_data, xs_rz, ys_rz, ro, zo, msing, psi_sing, q_sing = h5open(h5path, "r") do fid - read(fid["Equilibrium/Geometry/rcoords"]), read(fid["Equilibrium/Geometry/offset"]), - read(fid["Equilibrium/Geometry/psi"]), read(fid["Equilibrium/Geometry/theta"]), - read(fid["Equilibrium/R_axis"]), read(fid["Equilibrium/Z_axis"]), - read(fid["SingularSurfaces/rational_count"]), read(fid["SingularSurfaces/rational_psi"]), read(fid["SingularSurfaces/rational_q"]) - end +function plot_flux_surfaces(src::Source; n_psi=11, n_theta=18, save_path=nothing) + _require(plot_flux_surfaces, src) + p = _plot_flux_surfaces(_geometry_data(src); n_psi, n_theta) + isnothing(save_path) || savefig(p, save_path) + return p +end + +function _plot_flux_surfaces(d; n_psi=11, n_theta=18) + xs_rz, ys_rz, psi_sing, q_sing = d.psi, d.theta, d.sing_psi, d.sing_q n_psi_grid = length(xs_rz) n_theta_grid = length(ys_rz) - - # Build R, Z on the full nodal grid - R_grid = Matrix{Float64}(undef, n_psi_grid, n_theta_grid) - Z_grid = Matrix{Float64}(undef, n_psi_grid, n_theta_grid) - for ipsi in 1:n_psi_grid - for itheta in 1:n_theta_grid - rfac = sqrt(max(0.0, rcoords[ipsi, itheta])) - eta = 2π * (ys_rz[itheta] + offset_data[ipsi, itheta]) - R_grid[ipsi, itheta] = ro + rfac * cos(eta) - Z_grid[ipsi, itheta] = zo + rfac * sin(eta) - end - end + R_grid, Z_grid = _rz_grid(d) p = plot(; title="Flux surfaces in (R, Z)", @@ -216,7 +275,7 @@ function plot_flux_surfaces(h5path; n_psi=11, n_theta=18, save_path=nothing) color=:tomato, linewidth=0.8, label=label) end # Draw rational surface flux contours in red, annotated at the bottom (min Z) of each - for s in 1:msing + for s in eachindex(psi_sing) idx = argmin(abs.(xs_rz .- psi_sing[s])) q_label = abs(q_sing[s] - round(q_sing[s])) < 0.05 ? "q=$(round(Int, q_sing[s]))" : "q=$(round(q_sing[s], digits=2))" @@ -227,8 +286,6 @@ function plot_flux_surfaces(h5path; n_psi=11, n_theta=18, save_path=nothing) annotate!(p, R_grid[idx, itheta_bot], Z_grid[idx, itheta_bot], text(" $q_label", 7, :left, :red)) end - - isnothing(save_path) || savefig(p, save_path) return p end @@ -237,8 +294,8 @@ end Plot Grad-Shafranov error vs ψ_N (log scale) for several θ slices, with the flux-surface-integrated error overplotted as a thick black line. Reads `gse.h5` -(and optionally `gsei.h5`) from the same directory as `h5path`. Returns `nothing` -if `gse.h5` is not found (requires `diagnose_src = true` in the equilibrium configuration). +(and optionally `gsei.h5`) from the same directory as `h5path`, written when the equilibrium +configuration sets `diagnose_src = true`. ### Arguments @@ -251,11 +308,11 @@ if `gse.h5` is not found (requires `diagnose_src = true` in the equilibrium conf ### Returns -A `Plots.jl` plot object, or `nothing` if `gse.h5` is absent. +A `Plots.jl` plot object. """ -function plot_gse_by_theta(h5path; n_theta_lines=8, save_path=nothing) +function plot_gse_by_theta(h5path::AbstractString; n_theta_lines=8, save_path=nothing) + _require(plot_gse_by_theta, h5path) gse_path = joinpath(dirname(h5path), "gse.h5") - isfile(gse_path) || return nothing gse_data = h5open(gse_path, "r") do fid read(fid["gse_data"]) # (npsi, ntheta, 7): cols = θ, ψ, flux1, flux2, source, total, error @@ -297,8 +354,8 @@ end plot_gse_integrated(h5path; save_path=nothing) Plot the flux-surface-integrated Grad-Shafranov error (log scale) vs ψ_N, reading from -`gsei.h5` in the same directory as `h5path`. Returns `nothing` if `gsei.h5` is not found -(requires `diagnose_src = true` in the equilibrium configuration). +`gsei.h5` in the same directory as `h5path`, written when the equilibrium configuration sets +`diagnose_src = true`. ### Arguments @@ -310,11 +367,11 @@ Plot the flux-surface-integrated Grad-Shafranov error (log scale) vs ψ_N, readi ### Returns -A `Plots.jl` plot object, or `nothing` if `gsei.h5` is absent. +A `Plots.jl` plot object. """ -function plot_gse_integrated(h5path; save_path=nothing) +function plot_gse_integrated(h5path::AbstractString; save_path=nothing) + _require(plot_gse_integrated, h5path) gsei_path = joinpath(dirname(h5path), "gsei.h5") - isfile(gsei_path) || return nothing xs, errlogi = h5open(gsei_path, "r") do fid read(fid["xs"]), read(fid["errlogi"]) @@ -335,7 +392,7 @@ function plot_gse_integrated(h5path; save_path=nothing) end """ - plot_equilibrium_summary(h5path; save_path=nothing) + plot_equilibrium_summary(src; save_path=nothing) Summary of equilibrium profiles and geometry. The (R, Z) flux surface plot occupies the full left column. Profile plots are stacked in the right column: @@ -344,12 +401,12 @@ full left column. Profile plots are stacked in the right column: - μ₀p(ψ) pressure profile (`plot_pressure_profile`) - 2πF(ψ) toroidal field function (`plot_f_profile`) -If `gse.h5` is present (requires `diagnose_src = true`), a combined Grad-Shafranov error -panel (θ slices + integrated, log scale) is appended at the bottom of the right column. +When the source is a `gpec.h5` path with `gse.h5` beside it (requires `diagnose_src = true`), a combined +Grad-Shafranov error panel (θ slices + integrated, log scale) is appended at the bottom of the right column. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file + - `src`: A `PlasmaEquilibrium`, a `ForceFreeStatesResult`, or the path of a `gpec.h5` ### Keyword arguments @@ -359,20 +416,18 @@ panel (θ slices + integrated, log scale) is appended at the bottom of the right A `Plots.jl` plot object. """ -function plot_equilibrium_summary(h5path; save_path=nothing) - q0, q95, betat, betan, kappa, li1 = h5open(h5path, "r") do fid - read(fid["Equilibrium/q_axis"]), read(fid["Equilibrium/q_95"]), - read(fid["Equilibrium/beta_t"]), read(fid["Equilibrium/beta_N"]), - read(fid["Equilibrium/kappa"]), read(fid["Equilibrium/l_i_1"]) - end +function plot_equilibrium_summary(src::Source; save_path=nothing) + _require(plot_equilibrium_summary, src) + (; q0, q95, betat, betan, kappa, li1) = _global_data(src) title_str = "q0=$(round(q0,digits=2)) q95=$(round(q95,digits=2)) βₜ=$(round(betat,digits=3)) βₙ=$(round(betan,digits=3)) κ=$(round(kappa,digits=2)) li1=$(round(li1,digits=3))" - p_rz = plot_flux_surfaces(h5path) - p_q = plot_qprofile(h5path; show_singular=true) - p_pres = plot_pressure_profile(h5path) - p_f = plot_f_profile(h5path) - p_gse = plot_gse_by_theta(h5path) # includes integrated overplot; nothing if absent + profiles = _profile_data(src) + p_rz = _plot_flux_surfaces(_geometry_data(src)) + p_q = _plot_qprofile(profiles; show_singular=true) + p_pres = _plot_pressure_profile(profiles) + p_f = _plot_f_profile(profiles) + p_gse = src isa AbstractString && _can(plot_gse_by_theta, capabilities(src)) ? plot_gse_by_theta(src) : nothing # includes integrated overplot # Suppress x-axis labels/ticks on all but the bottom profile panel — they share the # same ψ_N axis and labeling every panel wastes vertical space. @@ -397,4 +452,12 @@ function plot_equilibrium_summary(h5path; save_path=nothing) return p end +register!(plot_equilibrium_summary, (PlasmaEquilibrium, ForceFreeStatesResult, AbstractString), (:equilibrium,), "flux surfaces and profiles") +register!(plot_qprofile, (PlasmaEquilibrium, ForceFreeStatesResult, AbstractString), (:equilibrium,), "q against ψ_N with the rational surfaces") +register!(plot_pressure_profile, (PlasmaEquilibrium, ForceFreeStatesResult, AbstractString), (:equilibrium,), "μ₀p against ψ_N") +register!(plot_f_profile, (PlasmaEquilibrium, ForceFreeStatesResult, AbstractString), (:equilibrium,), "2πF against ψ_N") +register!(plot_flux_surfaces, (PlasmaEquilibrium, ForceFreeStatesResult, AbstractString), (:equilibrium,), "flux surfaces and θ spokes in (R, Z)") +register!(plot_gse_by_theta, (AbstractString,), (:gse_files,), "Grad-Shafranov error against ψ_N per θ slice") +register!(plot_gse_integrated, (AbstractString,), (:gse_files,), "flux-surface-integrated Grad-Shafranov error") + end # module Equilibrium diff --git a/src/Analysis/ErrorFields.jl b/src/Analysis/ErrorFields.jl index ea980b400..635cc4745 100644 --- a/src/Analysis/ErrorFields.jl +++ b/src/Analysis/ErrorFields.jl @@ -17,10 +17,11 @@ using TOML import ...ErrorFields as EF import ...PerturbedEquilibrium as PE import ...ForcingTerms as FT +using ..Analysis: register!, _require, _can, _missing, _entry, _why, capabilities """ A label paired with something to plot: a `gpec.h5` path, or a result already in memory -(`ErrorFields.SensitivityTable`, `CoilSensitivities`, `MonteCarloResult`, `RiskResult`, or a +(`ErrorFields.ErrorFieldsResult`, `SensitivityTable`, `CoilSensitivities`, `MonteCarloResult`, `RiskResult`, or a vector of `CoilOverlap`). Analysing coil geometry that was never part of a run produces the latter, so the plots accept both. """ @@ -29,7 +30,7 @@ const Sources = AbstractVector{<:Pair{String,<:Any}} """ A single unlabelled source: a `gpec.h5` path, or one in-memory result. """ -const SingleSource = Union{AbstractString,EF.SensitivityTable,EF.CoilSensitivities,EF.MonteCarloResult,EF.RiskResult,AbstractVector{EF.CoilOverlap}} +const SingleSource = Union{AbstractString,EF.ErrorFieldsResult,EF.SensitivityTable,EF.CoilSensitivities,EF.MonteCarloResult,EF.RiskResult,AbstractVector{EF.CoilOverlap}} _sources(source) = [_default_label(source) => source] _default_label(h5path::AbstractString) = basename(dirname(abspath(h5path))) @@ -78,6 +79,30 @@ rather than separately, which would leave each half unable to find the other. _load(pair::Tuple{EF.MonteCarloResult,EF.RiskResult}) = _load(merge(_present(_load(pair[1])), _present(_load(pair[2])))) +# The whole assessment: the dominant-mode table and mode labels as the writer records them, from the plasma response it was made on. +function _load(r::EF.ErrorFieldsResult) + resp = r.response + ffs = resp.ffs + mn_index = hcat([(i - 1) % ffs.mpert + ffs.mlow for i in 1:ffs.numpert_total], [(i - 1) ÷ ffs.mpert + ffs.nlow for i in 1:ffs.numpert_total]) + # The dominant-mode table's shift and tilt sensitivities replace the full per-mode arrays, as in the file. + parts = Any[_present(_load(r.sensitivities)), _present(_load(EF.sensitivity_table(r.sensitivities, r.dominant; mode=1)))] + r.monte_carlo === nothing || push!(parts, _present(_load(r.monte_carlo))) + r.risk === nothing || push!(parts, _present(_load(r.risk))) + sc = r.scan + sc === nothing || push!(parts, + (; scan_scale=sc.tolerance_scale, scan_locking_probability_percent=sc.locking_probability_percent, + scan_locking_probability_efc_percent=sc.locking_probability_efc_percent, scan_locking_probability_spread_percent=sc.locking_probability_spread_percent, + scan_locking_probability_efc_spread_percent=sc.locking_probability_efc_spread_percent)) + ctrl = r.settings.ctrl + push!( + parts, + _present((; tolerances=r.tolerances, fd_step_shift_m=ctrl.fd_step_shift_m, fd_step_tilt_deg=ctrl.fd_step_tilt_deg, mn_index, + dominant_v=isempty(resp.dominant_singular_values) ? nothing : resp.dominant_right_singular_vectors, + singular_values=isempty(resp.dominant_singular_values) ? nothing : resp.dominant_singular_values)) + ) + return _load(merge(parts...)) +end + _load(ovs::AbstractVector{EF.CoilOverlap}) = _load((; coil_names=[o.coil_name for o in ovs], delta_as_designed=[o.delta for o in ovs], resonant_fraction_percent=[o.resonant_fraction_percent for o in ovs], b_t0=isempty(ovs) ? nothing : first(ovs).b_t0)) @@ -144,6 +169,10 @@ _centers(edges) = (edges[1:end-1] .+ edges[2:end]) ./ 2 _step_series(m, a) = (vcat(m[1] - 1, m, m[end] + 1), vcat(0.0, a, 0.0)) _empty(msg) = plot(; title=msg, legend=false) +# A registered source (a path or the assessment) must hold the data a single-source plot draws; other in-memory sources are what they are. +_check(fn, src::Union{AbstractString,EF.ErrorFieldsResult}) = _require(fn, src) +_check(fn, src) = nothing + function _save(p, save_path) if save_path !== nothing Plots.savefig(p, save_path) @@ -256,7 +285,7 @@ function _sensitivity_values(lbl, d, names, quantity::Symbol) vals === nothing && throw(ArgumentError("source \"$lbl\" carries no $field, so quantity=:$quantity cannot be drawn from it")) return abs.(vals[idx]) end -plot_coil_sensitivities(source::SingleSource; kwargs...) = plot_coil_sensitivities(_sources(source); kwargs...) +plot_coil_sensitivities(source::SingleSource; kwargs...) = (_check(plot_coil_sensitivities, source); plot_coil_sensitivities(_sources(source); kwargs...)) """ plot_tolerance_pdf(sources; corrected=true, normalize=false, xscale=:identity, show_batches=false, save_path=nothing) @@ -295,7 +324,7 @@ function plot_tolerance_pdf(sources::Sources; corrected::Bool=true, normalize::B any_data || return _empty("No ErrorFields/MonteCarlo data — run with a tolerance_file") return _save(p, save_path) end -plot_tolerance_pdf(source::SingleSource; kwargs...) = plot_tolerance_pdf(_sources(source); kwargs...) +plot_tolerance_pdf(source::SingleSource; kwargs...) = (_check(plot_tolerance_pdf, source); plot_tolerance_pdf(_sources(source); kwargs...)) """ plot_overlap_phasors(sources; coils=nothing, save_path=nothing) @@ -347,7 +376,7 @@ function plot_overlap_phasors(sources::Sources; coils=nothing, save_path=nothing length(panels) == 1 && return _save(panels[1], save_path) return _save(plot(panels...; layout=(1, length(panels)), size=(900 * length(panels), 560)), save_path) end -plot_overlap_phasors(source::SingleSource; kwargs...) = plot_overlap_phasors(_sources(source); kwargs...) +plot_overlap_phasors(source::SingleSource; kwargs...) = (_check(plot_overlap_phasors, source); plot_overlap_phasors(_sources(source); kwargs...)) """ plot_tolerance_budget(h5path; psi_low=0.0, psi_high=CORE_PSI_HIGH, mode=1, tolerance_scale=1.0, sort=:total, top=nothing, save_path=nothing) @@ -472,7 +501,7 @@ function plot_linearity_residuals(sources::Sources; at::Symbol=:tolerance, coils length(panels) == 1 && return _save(panels[1], save_path) return _save(plot(panels...; layout=(1, length(panels)), size=(900 * length(panels), 520)), save_path) end -plot_linearity_residuals(source::SingleSource; kwargs...) = plot_linearity_residuals(_sources(source); kwargs...) +plot_linearity_residuals(source::SingleSource; kwargs...) = (_check(plot_linearity_residuals, source); plot_linearity_residuals(_sources(source); kwargs...)) """ plot_linearity_check(check::EF.LinearityCheck; save_path=nothing) @@ -593,7 +622,7 @@ function plot_locking_risk(sources::Sources; corrected::Bool=true, target_percen target_percent === nothing || hline!(p, [target_percent]; ls=:dash, c=:gray, label="target $(target_percent) %") return _save(p, save_path) end -plot_locking_risk(source::SingleSource; kwargs...) = plot_locking_risk(_sources(source); kwargs...) +plot_locking_risk(source::SingleSource; kwargs...) = (_check(plot_locking_risk, source); plot_locking_risk(_sources(source); kwargs...)) """ plot_threshold_scaling(sources; save_path=nothing) @@ -623,7 +652,7 @@ function plot_threshold_scaling(sources::Sources; save_path=nothing) any_data || return _empty("No ErrorFields/Risk data — run with an [ErrorFields.scenario] table") return _save(p, save_path) end -plot_threshold_scaling(source::SingleSource; kwargs...) = plot_threshold_scaling(_sources(source); kwargs...) +plot_threshold_scaling(source::SingleSource; kwargs...) = (_check(plot_threshold_scaling, source); plot_threshold_scaling(_sources(source); kwargs...)) """ plot_dominant_mode_spectrum(sources; mode=1, save_path=nothing) @@ -651,7 +680,7 @@ function plot_dominant_mode_spectrum(sources::Sources; mode::Int=1, save_path=no any_data || return _empty("No PerturbedEquilibrium/SingularCoupling/DominantMode data") return _save(p, save_path) end -plot_dominant_mode_spectrum(source::SingleSource; kwargs...) = plot_dominant_mode_spectrum(_sources(source); kwargs...) +plot_dominant_mode_spectrum(source::SingleSource; kwargs...) = (_check(plot_dominant_mode_spectrum, source); plot_dominant_mode_spectrum(_sources(source); kwargs...)) """ plot_applied_spectra(ctx, overlaps; normalize=true, save_path=nothing) @@ -987,17 +1016,35 @@ function plot_needed_current(mc::EF.MonteCarloResult, array::EF.CoilOverlap; cou end """ - plot_error_field_summary(h5path; save_path=nothing) + plot_error_field_summary(src; save_path=nothing) -Four panels of a run's error-field assessment: coil sensitivities to shift, the tolerance -Monte Carlo distributions, the overlap distribution against the threshold, and the locking -risk against tolerance scale. Panels whose data the run did not produce say so. +Panels of a run's error-field assessment, as many as it produced: coil sensitivities to shift, the tolerance +Monte Carlo distributions, the overlap distribution against the threshold, and the locking risk against +tolerance scale. `src` is an `ErrorFields.ErrorFieldsResult` or the path of a `gpec.h5` it was saved to. """ -function plot_error_field_summary(h5path::AbstractString; save_path=nothing) - src = _sources(h5path) - p = plot(plot_coil_sensitivities(src), plot_tolerance_pdf(src), plot_threshold_scaling(src), plot_locking_risk(src); - layout=(2, 2), size=(1400, 1000)) +function plot_error_field_summary(source::Union{AbstractString,EF.ErrorFieldsResult}; save_path=nothing) + caps = capabilities(source) + src = _sources(source) + panels = (plot_coil_sensitivities, plot_tolerance_pdf, plot_threshold_scaling, plot_locking_risk) + drawn = [fn(src) for fn in panels if _can(fn, caps)] + if isempty(drawn) + missing = unique(reduce(vcat, (_missing(_entry(fn), caps) for fn in panels))) + throw(ArgumentError("plot_error_field_summary: no panel can be drawn: $(_why(missing))")) + end + n = length(drawn) + n == 3 && push!(drawn, plot(; framestyle=:none)) + p = n <= 2 ? plot(drawn...; layout=(1, n), size=(700n, 500)) : plot(drawn...; layout=(2, 2), size=(1400, 1000)) return _save(p, save_path) end +const _EF_SOURCES = (EF.ErrorFieldsResult, AbstractString) +register!(plot_error_field_summary, _EF_SOURCES, (), "sensitivities, tolerance Monte Carlo, threshold and locking risk") +register!(plot_coil_sensitivities, _EF_SOURCES, (:ef_sensitivities,), "per-coil dominant-mode sensitivities") +register!(plot_overlap_phasors, _EF_SOURCES, (:ef_sensitivities,), "as-designed overlap δ of each coil set in the complex plane") +register!(plot_linearity_residuals, _EF_SOURCES, (:ef_sensitivities,), "linearity residuals of the coil sensitivities") +register!(plot_dominant_mode_spectrum, _EF_SOURCES, (:dominant_mode,), "spectrum of the dominant resonant-coupling mode") +register!(plot_tolerance_pdf, _EF_SOURCES, (:ef_monte_carlo,), "Monte Carlo distribution of |δ|") +register!(plot_threshold_scaling, _EF_SOURCES, (:ef_monte_carlo, :ef_risk), "|δ| distribution against the penetration threshold") +register!(plot_locking_risk, _EF_SOURCES, (:ef_scan,), "locking risk against tolerance scale") + end # module ErrorFields diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index 47e839244..b36636e25 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -1,8 +1,9 @@ """ ForceFreeStates -Post-processing and visualization functions for ForceFreeStates (DCON-style ideal MHD stability) -results stored in GPEC HDF5 output files. +Post-processing and visualization functions for ForceFreeStates (DCON-style ideal MHD stability) results. Each plot +takes the `ForceFreeStatesResult` in memory or the path of a `gpec.h5` it was saved to, and throws an `ArgumentError` +naming what to change in the run when the source lacks the data it draws (see `Analysis.available_plots`). """ module ForceFreeStates @@ -12,15 +13,45 @@ using Plots using Printf using Statistics: quantile +using ..Analysis: register!, _require, _can, _missing, _entry, _why, capabilities +import ...ForceFreeStates: ForceFreeStatesResult +import ...LocalStability: LocalStabilityResult + +const Source = Union{AbstractString,ForceFreeStatesResult} + +const _FWD = "ForceFreeStates/Solutions/ForwardIntegration" +const _GAL = "ForceFreeStates/Solutions/GalerkinIntegration" +const _FBS = "ForceFreeStates/FreeBoundaryStability" + +_label(h5path::AbstractString) = h5path +_label(::ForceFreeStatesResult) = "in memory" + +# Least-stable total energy, or `nothing` without free-boundary energies. +_first_energy(et) = isempty(et) ? nothing : et[1] + +function _mode_displacement_data(r::ForceFreeStatesResult) + sol = r.solution + return (; mlow=r.mlow, psi=sol.psi_store, xi=sol.u_store[:, 1, 1, :], dW=r.free_boundary === nothing ? nothing : _first_energy(r.free_boundary.et)) +end + +function _mode_displacement_data(h5path::AbstractString) + h5open(h5path, "r") do fid + # Forward integration stores its profiles in ForwardIntegration, a matched Galerkin solve in GalerkinIntegration. + group = isempty(fid["$_FWD/xi_psi"]) ? _GAL : _FWD + (; mlow=read(fid["Info/mlow"]), psi=read(fid["$group/psi"]), xi=read(fid["$group/xi_psi"])[:, 1, :], + dW=_first_energy(read(fid["$_FBS/eigenmode_energies"]))) + end +end + """ - plot_mode_displacement(h5path; modes=1:5, save_path=nothing) + plot_mode_displacement(src; modes=1:5, save_path=nothing) Plot |ξ_ψ| vs ψ_N for the least stable eigenmode, showing one curve per requested poloidal mode number m. The title includes the first eigenvalue dW = et[1]. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file (e.g. `"gpec.h5"`) + - `src`: A `ForceFreeStatesResult` with stored ξ profiles, or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -31,17 +62,18 @@ poloidal mode number m. The title includes the first eigenvalue dW = et[1]. A `Plots.jl` plot object. """ -function plot_mode_displacement(h5path; modes=1:5, save_path=nothing) - mlow, xi_psi, psi, et = h5open(h5path, "r") do fid - read(fid["Info/mlow"]), read(fid["ForceFreeStates/Solutions/ForwardIntegration/xi_psi"]), - read(fid["ForceFreeStates/Solutions/ForwardIntegration/psi"]), read(fid["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"]) - end +function plot_mode_displacement(src::Source; modes=1:5, save_path=nothing) + _require(plot_mode_displacement, src) + p = _plot_mode_displacement(_mode_displacement_data(src); modes) + isnothing(save_path) || savefig(p, save_path) + return p +end - mpert = size(xi_psi, 1) - mhigh = mlow + mpert - 1 - dW = isempty(et) ? nothing : et[1] - title_str = isnothing(dW) ? "Least stable mode" : - "Least stable mode, δW = $(round(real(dW), sigdigits=4))" +function _plot_mode_displacement(d; modes=1:5) + mpert = size(d.xi, 1) + mhigh = d.mlow + mpert - 1 + title_str = isnothing(d.dW) ? "Least stable mode" : + "Least stable mode, δW = $(round(real(d.dW), sigdigits=4))" p = plot(; xlims=(0, 1), @@ -52,24 +84,24 @@ function plot_mode_displacement(h5path; modes=1:5, save_path=nothing) bottom_margin=5Plots.mm ) for m in modes - mlow <= m <= mhigh || continue - plot!(p, psi, abs.(xi_psi[m-mlow+1, 1, :]); label="m=$m") + d.mlow <= m <= mhigh || continue + plot!(p, d.psi, abs.(d.xi[m-d.mlow+1, :]); label="m=$m") end - - isnothing(save_path) || savefig(p, save_path) return p end +_criterion_data(r::ForceFreeStatesResult) = (; psi=r.diagnostics.psi_store, crit=r.diagnostics.crit_store) +_criterion_data(h5path::AbstractString) = h5open(fid -> (; psi=read(fid["$_FWD/psi"]), crit=read(fid["$_FWD/crit"])), h5path, "r") """ - plot_fixed_boundary_stability_criterion(h5path; save_path=nothing) + plot_fixed_boundary_stability_criterion(src; save_path=nothing) -Plot the stability criterion (smallest eigenvalue of W⁻¹, `crit`) vs ψ_N. +Plot the critical determinant (smallest eigenvalue of W⁻¹, `crit`) vs ψ_N. A sign change in `crit` during integration indicates an ideal fixed-boundary instability. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file + - `src`: A Forward or Riccati `ForceFreeStatesResult`, or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -79,13 +111,16 @@ A sign change in `crit` during integration indicates an ideal fixed-boundary ins A `Plots.jl` plot object. """ -function plot_fixed_boundary_stability_criterion(h5path; save_path=nothing) - psi, crit = h5open(h5path, "r") do fid - read(fid["ForceFreeStates/Solutions/ForwardIntegration/psi"]), read(fid["ForceFreeStates/Solutions/ForwardIntegration/crit"]) - end +function plot_fixed_boundary_stability_criterion(src::Source; save_path=nothing) + _require(plot_fixed_boundary_stability_criterion, src) + p = _plot_fixed_boundary_stability_criterion(_criterion_data(src)) + isnothing(save_path) || savefig(p, save_path) + return p +end +function _plot_fixed_boundary_stability_criterion(d) p = plot( - psi, crit; + d.psi, d.crit; xlims=(0, 1), xlabel="Norm. Poloidal Flux", ylabel="|Dᶜ|", @@ -95,25 +130,26 @@ function plot_fixed_boundary_stability_criterion(h5path; save_path=nothing) bottom_margin=5Plots.mm ) hline!(p, [0.0]; linestyle=:dash, color=:black, label=nothing) - - isnothing(save_path) || savefig(p, save_path) return p end +_energy_eigenvector_data(r::ForceFreeStatesResult) = (; wt=r.free_boundary.wt, psio=r.equil.psio, mlow=r.mlow) +_energy_eigenvector_data(h5path::AbstractString) = + h5open(fid -> (; wt=read(fid["$_FBS/W_freeboundary_eigenmodes"]), psio=read(fid["Equilibrium/psi_total"]), mlow=read(fid["Info/mlow"])), h5path, "r") + """ - plot_energy_eigenvectors(h5path; matrix_type=:total, save_path=nothing) + plot_energy_eigenvectors(src; matrix_type=:total, save_path=nothing) Heatmap of energy eigenvector magnitudes vs (m, mode index). -Only `matrix_type=:total` is supported (the total energy eigenvector matrix `Wₜ` is stored in -`FreeBoundaryStability/W_freeboundary_eigenmodes`). Plasma and vacuum eigenvectors are not -stored separately in the HDF5 output. +Only `matrix_type=:total` is supported (the total energy eigenvector matrix `Wₜ`). Plasma and vacuum eigenvectors +are not stored separately. Eigenvectors are scaled by χ₁ = 2π ψ₀ × 10⁻³ to match GPEC conventions. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file with vacuum data (`vac_flag = true`) + - `src`: A `ForceFreeStatesResult` with free-boundary energies (`vac_flag = true`), or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -124,24 +160,24 @@ Eigenvectors are scaled by χ₁ = 2π ψ₀ × 10⁻³ to match GPEC convention A `Plots.jl` plot object. """ -function plot_energy_eigenvectors(h5path; matrix_type=:total, save_path=nothing) +function plot_energy_eigenvectors(src::Source; matrix_type=:total, save_path=nothing) matrix_type == :total || - error("matrix_type=$matrix_type not supported; only :total has eigenvector matrix stored in HDF5 (ep/ev are eigenvalue vectors, not matrices)") - - wt, psio, mlow = h5open(h5path, "r") do fid - read(fid["ForceFreeStates/FreeBoundaryStability/W_freeboundary_eigenmodes"]), read(fid["Equilibrium/psi_total"]), read(fid["Info/mlow"]) - end - - isempty(wt) && error("No vacuum data in $h5path; rerun with vac_flag = true") + throw(ArgumentError("matrix_type=$matrix_type not supported; only :total has an eigenvector matrix (ep/ev are eigenvalue vectors, not matrices)")) + _require(plot_energy_eigenvectors, src) + p = _plot_energy_eigenvectors(_energy_eigenvector_data(src)) + isnothing(save_path) || savefig(p, save_path) + return p +end - chi1 = 2π * psio - wt_scaled = wt * (chi1 * 1e-3) +function _plot_energy_eigenvectors(d) + chi1 = 2π * d.psio + wt_scaled = d.wt * (chi1 * 1e-3) nmn = size(wt_scaled, 1) nmodes = size(wt_scaled, 2) - m_vals = (0:(nmn-1)) .+ mlow + m_vals = (0:(nmn-1)) .+ d.mlow - p = heatmap( + return heatmap( m_vals, 1:nmodes, abs.(wt_scaled'); xlabel="Poloidal Harmonic", ylabel="Eigenmode Index", @@ -151,13 +187,23 @@ function plot_energy_eigenvectors(h5path; matrix_type=:total, save_path=nothing) right_margin=20Plots.mm, bottom_margin=5Plots.mm ) +end - isnothing(save_path) || savefig(p, save_path) - return p +function _edge_scan_data(r::ForceFreeStatesResult) + es = r.diagnostics.edge_scan + return (; q=es.q, et=es.total_eigenvalue, ep=es.plasma_energy, ev=es.vacuum_energy, evonly=es.vacuum_eigenvalue, qlim=r.qlim, label=_label(r)) +end + +function _edge_scan_data(h5path::AbstractString) + h5open(h5path, "r") do fid + (; q=read(fid["ForceFreeStates/EdgeScan/q"]), et=read(fid["ForceFreeStates/EdgeScan/total_energy"]), + ep=read(fid["ForceFreeStates/EdgeScan/plasma_energy"]), ev=read(fid["ForceFreeStates/EdgeScan/vacuum_energy"]), + evonly=read(fid["ForceFreeStates/EdgeScan/vacuum_eigenvalue"]), qlim=read(fid["Info/qlim"]), label=_label(h5path)) + end end """ - plot_edge_stability_scan(h5path; save_path=nothing, ylims=(-2, 3), kwargs...) + plot_edge_stability_scan(src; save_path=nothing, ylims=(-2, 3), kwargs...) Plot the edge stability scan energy components (et, ep, ev, evonly) vs ψ_N. @@ -177,7 +223,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 `dW_edge_scan_start < psilim` + - `src`: A `ForceFreeStatesResult` produced with `dW_edge_scan_start < psilim`, or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -187,55 +233,41 @@ A horizontal dashed line at zero marks the stability boundary. A vertical dashed ### Returns -A `Plots.jl` plot object, or `nothing` if no `EdgeScan/` group is present in the file. +A `Plots.jl` plot object. """ -function plot_edge_stability_scan(h5path; save_path=nothing, ylims=(-2, 3), kwargs...) - has_scan, q, et, ep, ev, evonly, qlim = h5open(h5path, "r") do fid - if !haskey(fid, "ForceFreeStates/EdgeScan/psi") - return false, Float64[], ComplexF64[], ComplexF64[], ComplexF64[], Float64[], NaN - end - true, - read(fid["ForceFreeStates/EdgeScan/q"]), - read(fid["ForceFreeStates/EdgeScan/total_energy"]), - read(fid["ForceFreeStates/EdgeScan/plasma_energy"]), - read(fid["ForceFreeStates/EdgeScan/vacuum_energy"]), - read(fid["ForceFreeStates/EdgeScan/vacuum_eigenvalue"]), - read(fid["Info/qlim"]) - end - - if !has_scan - @warn "No edge_scan group in $h5path. Run with dW_edge_scan_start < psilim to generate it." - return nothing - end +function plot_edge_stability_scan(src::Source; save_path=nothing, ylims=(-2, 3), kwargs...) + _require(plot_edge_stability_scan, src) + d = _edge_scan_data(src) + q, qlim = d.q, d.qlim kw_re = (xlabel="q", label="Re", ylims=ylims, kwargs...) kw_im = (xlabel="q", label="Im", ls=:dash, ylims=ylims, kwargs...) vl_kw = (color=:gray, lw=1, ls=:dot, label=false) hl_kw = (color=:black, lw=1, ls=:dash, label=false) - p_et = plot(q, real.(et); ylabel="Total Energy", title="Edge Stability Scan: δW vs q", kw_re...) - plot!(p_et, q, imag.(et); kw_im...) + p_et = plot(q, real.(d.et); ylabel="Total Energy", title="Edge Stability Scan: δW vs q", kw_re...) + plot!(p_et, q, imag.(d.et); kw_im...) hline!(p_et, [0.0]; hl_kw...) vline!(p_et, [qlim]; vl_kw...) - p_ep = plot(q, real.(ep); ylabel="Plasma Energy", kw_re...) - plot!(p_ep, q, imag.(ep); kw_im...) + p_ep = plot(q, real.(d.ep); ylabel="Plasma Energy", kw_re...) + plot!(p_ep, q, imag.(d.ep); kw_im...) hline!(p_ep, [0.0]; hl_kw...) vline!(p_ep, [qlim]; vl_kw...) - p_ev = plot(q, real.(ev); ylabel="Vacuum Energy", kw_re...) - plot!(p_ev, q, imag.(ev); kw_im...) + p_ev = plot(q, real.(d.ev); ylabel="Vacuum Energy", kw_re...) + plot!(p_ev, q, imag.(d.ev); kw_im...) hline!(p_ev, [0.0]; hl_kw...) vline!(p_ev, [qlim]; vl_kw...) - p_evonly = plot(q, evonly; ylabel="Min Vac. Eigenvalue", legend=false, xlabel="q", ylims=ylims, kwargs...) + p_evonly = plot(q, d.evonly; ylabel="Min Vac. Eigenvalue", legend=false, xlabel="q", ylims=ylims, kwargs...) hline!(p_evonly, [0.0]; hl_kw...) vline!(p_evonly, [qlim]; vl_kw...) p = plot(p_et, p_ep, p_ev, p_evonly; layout=(4, 1), size=(900, 1100), - plot_title="Edge stability scan: $h5path", + plot_title="Edge stability scan: $(d.label)", left_margin=12Plots.mm, bottom_margin=4Plots.mm) @@ -243,15 +275,21 @@ function plot_edge_stability_scan(h5path; save_path=nothing, ylims=(-2, 3), kwar return p end +const _ENERGY_FIELDS = Dict(:total => :et, :plasma => :ep, :vacuum => :ev) +const _ENERGY_DATASETS = Dict(:total => "eigenmode_energies", :plasma => "eigenmode_plasma_energies", :vacuum => "eigenmode_vacuum_energies") + +_eigenvalue_data(r::ForceFreeStatesResult, matrix_type) = getfield(r.free_boundary, _ENERGY_FIELDS[matrix_type]) +_eigenvalue_data(h5path::AbstractString, matrix_type) = h5open(fid -> read(fid["$_FBS/$(_ENERGY_DATASETS[matrix_type])"]), h5path, "r") + """ - plot_eigenvalues(h5path; matrix_type=:total, save_path=nothing) + plot_eigenvalues(src; matrix_type=:total, save_path=nothing) -Scatter plot of energy eigenvalues vs mode index. Points are colored red (unstable, Re > 0) -or green (stable, Re < 0), with a dashed reference line at zero. +Scatter plot of energy eigenvalues vs mode index. Points are colored red (unstable, Re < 0) +or blue (stable, Re > 0), with a dashed reference line at zero. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file with vacuum data (`vac_flag = true`) + - `src`: A `ForceFreeStatesResult` with free-boundary energies (`vac_flag = true`), or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -262,19 +300,10 @@ or green (stable, Re < 0), with a dashed reference line at zero. A `Plots.jl` plot object. """ -function plot_eigenvalues(h5path; matrix_type=:total, save_path=nothing) - dataset = Dict( - :total => "ForceFreeStates/FreeBoundaryStability/eigenmode_energies", - :plasma => "ForceFreeStates/FreeBoundaryStability/eigenmode_plasma_energies", - :vacuum => "ForceFreeStates/FreeBoundaryStability/eigenmode_vacuum_energies" - ) - haskey(dataset, matrix_type) || error("matrix_type must be :total, :plasma, or :vacuum") - - et = h5open(h5path, "r") do fid - read(fid[dataset[matrix_type]]) - end - - isempty(et) && error("No vacuum data in $h5path; rerun with vac_flag = true") +function plot_eigenvalues(src::Source; matrix_type=:total, save_path=nothing) + haskey(_ENERGY_FIELDS, matrix_type) || throw(ArgumentError("matrix_type must be :total, :plasma, or :vacuum")) + _require(plot_eigenvalues, src) + et = _eigenvalue_data(src, matrix_type) nmodes = length(et) ev_real = real.(et) @@ -297,19 +326,33 @@ function plot_eigenvalues(h5path; matrix_type=:total, save_path=nothing) return p end -""" - plot_delta_prime(h5path; save_path=nothing) +# Δ′_ii on the surfaces the Δ′ matrix spans: Galerkin solves across its in-band subset of the rational surfaces. +function _delta_prime_data(r::ForceFreeStatesResult) + dp = r.delta_prime.matrix + k = minimum(size(dp)) + g = r.galerkin + psi, q = g !== nothing && length(g.sing_psi) == k ? (g.sing_psi, g.sing_q) : ([s.psifac for s in r.surfaces[1:k]], [s.q for s in r.surfaces[1:k]]) + return (; psi, q, delta_prime=[dp[i, i] for i in 1:k]) +end + +function _delta_prime_data(h5path::AbstractString) + h5open(h5path, "r") do fid + dp = read(fid["SingularSurfaces/Delta_prime_matrix"]) + k = minimum(size(dp)) + group = haskey(fid, "$_GAL/rational_psi") && length(fid["$_GAL/rational_psi"]) == k ? _GAL : "SingularSurfaces" + (; psi=read(fid["$group/rational_psi"])[1:k], q=read(fid["$group/rational_q"])[1:k], delta_prime=[dp[i, i] for i in 1:k]) + end +end -Scatter plot of `Re(Δ')` per singular surface vs ψ_N, computed from the stored asymptotic -coefficients `ca_left` and `ca_right`. Points are colored red (tearing unstable, Re(Δ') > 0) -or blue (tearing stable). Integer-valued q rational surfaces are annotated. +""" + plot_delta_prime(src; save_path=nothing) -Δ' is computed as `(ca_right[resnum,resnum,2,s] - ca_left[resnum,resnum,2,s]) / (4π² ψ₀)`, -where `resnum` is the linear mode index of the (m,n) resonant pair at surface `s`. +Scatter plot of the diagonal `Re(Δ′_ii)` of the inter-surface Δ′ matrix (PEST-3 convention, ψ_N-referenced) vs ψ_N +of each rational surface it spans. Points are colored red (tearing unstable, Re(Δ′) > 0) or blue (tearing stable). ### Arguments - - `h5path`: Path to a GPEC HDF5 output file + - `src`: A Riccati or Galerkin `ForceFreeStatesResult`, or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -319,43 +362,22 @@ where `resnum` is the linear mode index of the (m,n) resonant pair at surface `s A `Plots.jl` plot object. """ -function plot_delta_prime(h5path; save_path=nothing) - msing, psi_sing, q_sing, ca_l, ca_r, psio, mn_index = h5open(h5path, "r") do fid - read(fid["SingularSurfaces/rational_count"]), read(fid["SingularSurfaces/rational_psi"]), read(fid["SingularSurfaces/rational_q"]), - read(fid["SingularSurfaces/ca_left"]), read(fid["SingularSurfaces/ca_right"]), - read(fid["Equilibrium/psi_total"]), read(fid["Info/mn_index"]) - end - - msing == 0 && return plot(; title="No singular surfaces found", legend=false) - # ca_left/ca_right are zero-extent sentinels on kinetic/galerkin-matched runs (never computed there). - isempty(ca_l) && return plot(; title="No asymptotic coefficients — ca_left/ca_right not computed for this run", legend=false) - - numpert_total = size(ca_l, 1) - chi1 = 2π * psio - - dp_vals = ComplexF64[] - for s in 1:msing - q_s = q_sing[s] - resnum = findfirst(1:numpert_total) do j - n_j = mn_index[j, 2] - n_j != 0 && abs(mn_index[j, 1] / n_j - q_s) < 1e-6 - end - if isnothing(resnum) - push!(dp_vals, 0.0 + 0.0im) - else - dp = (ca_r[resnum, resnum, 2, s] - ca_l[resnum, resnum, 2, s]) / (2π * chi1) - push!(dp_vals, dp) - end - end +function plot_delta_prime(src::Source; save_path=nothing) + _require(plot_delta_prime, src) + p = _plot_delta_prime(_delta_prime_data(src)) + isnothing(save_path) || savefig(p, save_path) + return p +end - dp_real = real.(dp_vals) +function _plot_delta_prime(d) + dp_real = real.(d.delta_prime) colors = [v > 0 ? :red : :steelblue for v in dp_real] p = scatter( - psi_sing, dp_real; + d.psi, dp_real; xlims=(0, 1), xlabel="Norm. Poloidal Flux", - ylabel="Re(Δ')", + ylabel="Re(Δ'_ii)", title="Tearing stability Δ'", legend=false, color=colors, @@ -365,23 +387,24 @@ function plot_delta_prime(h5path; save_path=nothing) bottom_margin=5Plots.mm ) hline!(p, [0.0]; linestyle=:dash, color=:black, label=nothing) - - isnothing(save_path) || savefig(p, save_path) return p end +_ballooning_data(r::LocalStabilityResult) = (; psi=r.ballooning_psi, alpha=r.alpha, alpha_crit=r.alpha_critical) +_ballooning_data(h5path::AbstractString) = + h5open(fid -> (; psi=read(fid["LocalStability/ballooning_psi"]), alpha=read(fid["LocalStability/alpha"]), alpha_crit=read(fid["LocalStability/alpha_critical"])), h5path, "r") + """ - plot_ballooning_alpha_boundary(h5path; save_path=nothing, psi_min=0.0) + plot_ballooning_alpha_boundary(src; save_path=nothing, psi_min=0.0) Plot the BALOO-style infinite-n ballooning stability diagram: the experimental pressure gradient α (solid) and the first stability boundary α_crit (dashed) versus normalized poloidal flux ψ_N. Surfaces where the experimental α lies above the boundary -are ballooning-unstable. Reads `LocalStability/ballooning_psi`, `LocalStability/alpha`, and -`LocalStability/alpha_critical` (populated when the deck has a `[LocalStability]` section). +are ballooning-unstable. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file (e.g. `"gpec.h5"`) + - `src`: A `LocalStabilityResult`, or the path of a `gpec.h5` written with a `[LocalStability]` section ### Keyword arguments @@ -393,13 +416,9 @@ are ballooning-unstable. Reads `LocalStability/ballooning_psi`, `LocalStability/ A `Plots.jl` plot object. """ -function plot_ballooning_alpha_boundary(h5path; save_path=nothing, psi_min=0.0) - psi, alpha, alpha_crit = h5open(h5path, "r") do fid - haskey(fid, "LocalStability/alpha") || return (Float64[], Float64[], Float64[]) - read(fid["LocalStability/ballooning_psi"]), read(fid["LocalStability/alpha"]), read(fid["LocalStability/alpha_critical"]) - end - - isempty(alpha) && return plot(; title="No local stability data (add a [LocalStability] section)", legend=false) +function plot_ballooning_alpha_boundary(src::Union{AbstractString,LocalStabilityResult}; save_path=nothing, psi_min=0.0) + _require(plot_ballooning_alpha_boundary, src) + d = _ballooning_data(src) p = plot(; xlims=(psi_min, 1), @@ -414,18 +433,31 @@ function plot_ballooning_alpha_boundary(h5path; save_path=nothing, psi_min=0.0) # NaN entries leave natural gaps over always-stable surfaces; isfinite masking would # wrongly connect across such gaps with a straight line. - plot!(p, psi, alpha; lw=2, color=:black, label="Experimental gradient") - plot!(p, psi, alpha_crit; lw=2, linestyle=:dash, color=:red, label="1st stability boundary") + plot!(p, d.psi, d.alpha; lw=2, color=:black, label="Experimental gradient") + plot!(p, d.psi, d.alpha_crit; lw=2, linestyle=:dash, color=:red, label="1st stability boundary") isnothing(save_path) || savefig(p, save_path) return p end +function _cond_fbar_data(r::ForceFreeStatesResult) + k = r.kinetic + return (; scan_psi=k.scan_psi, scan_cond=k.scan_cond, thr=k.scan_threshold, k_psi=[s.psifac for s in k.kinsing], + i_psi=[s.psifac for s in r.surfaces], i_q=[s.q for s in r.surfaces], kmsing=k.kmsing) +end + +function _cond_fbar_data(h5path::AbstractString) + h5open(h5path, "r") do fid + kg = fid["SingularSurfaces/Kinetic"] + (; scan_psi=read(kg["scan_psi"]), scan_cond=read(kg["scan_cond"]), thr=read(kg["scan_threshold"]), k_psi=read(kg["rational_psi"]), + i_psi=read(fid["SingularSurfaces/rational_psi"]), i_q=read(fid["SingularSurfaces/rational_q"]), kmsing=read(kg["rational_count"])) + end +end + """ - plot_cond_fbar(h5path; save_path=nothing, zoom=false) + plot_cond_fbar(src; save_path=nothing, zoom=false) -Plot `cond(F̄)` vs ψ from the kinetic-singular-surface scan stored in -`SingularSurfaces/Kinetic/` (populated when ForceFreeStates runs with +Plot `cond(F̄)` vs ψ from the kinetic-singular-surface scan (populated when ForceFreeStates runs with a `kinetic` source and `singfac_min > 0`). `F̄` is the kinetic Euler-Lagrange matrix formed by Schur-reducing the six @@ -443,7 +475,7 @@ kinetic resonances land relative to the ideal ones. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file produced with kinetic mode enabled + - `src`: A kinetic `ForceFreeStatesResult`, or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -454,27 +486,11 @@ kinetic resonances land relative to the ideal ones. ### Returns -A `Plots.jl` plot object, or `nothing` if no kinetic scan is stored in the file. +A `Plots.jl` plot object. """ -function plot_cond_fbar(h5path; save_path=nothing, zoom=false) - scan_psi, scan_cond, thr, k_psi, i_psi, i_q, kmsing = h5open(h5path, "r") do fid - if !(haskey(fid, "SingularSurfaces") && haskey(fid["SingularSurfaces"], "Kinetic")) - return Float64[], Float64[], 0.0, Float64[], Float64[], Float64[], 0 - end - kg = fid["SingularSurfaces/Kinetic"] - (read(kg["scan_psi"]), - read(kg["scan_cond"]), - read(kg["scan_threshold"]), - read(kg["rational_psi"]), - read(fid["SingularSurfaces/rational_psi"]), - read(fid["SingularSurfaces/rational_q"]), - read(kg["rational_count"])) - end - - if isempty(scan_psi) - @warn "No kinetic-singular-surface scan in $h5path — rerun with a kinetic source and singfac_min>0" - return nothing - end +function plot_cond_fbar(src::Source; save_path=nothing, zoom=false) + _require(plot_cond_fbar, src) + (; scan_psi, scan_cond, thr, k_psi, i_psi, i_q, kmsing) = _cond_fbar_data(src) finite_cond = filter(isfinite, scan_cond) data_max = isempty(finite_cond) ? 1.0 : maximum(finite_cond) @@ -527,22 +543,29 @@ function plot_cond_fbar(h5path; save_path=nothing, zoom=false) return p end +# The summary's panels in drawing order, each with the plot whose capabilities it needs. +_ffs_panels(src) = ( + plot_energy_eigenvectors => () -> _plot_energy_eigenvectors(_energy_eigenvector_data(src)), + plot_fixed_boundary_stability_criterion => () -> _plot_fixed_boundary_stability_criterion(_criterion_data(src)), + plot_mode_displacement => () -> _plot_mode_displacement(_mode_displacement_data(src)), + plot_delta_prime => () -> _plot_delta_prime(_delta_prime_data(src)) +) + """ - plot_ffs_summary(h5path; save_path=nothing) + plot_ffs_summary(src; save_path=nothing) -Four-panel summary of ForceFreeStates (DCON-style) stability results, combining: +Summary of ForceFreeStates (DCON-style) stability results, with each panel the source can supply: - Energy eigenvector heatmap (`plot_energy_eigenvectors`) - - Fixed-boundary stability criterion |D_c| vs ψ_N (`plot_stability_criterion`) - - Eigenvalue spectrum (`plot_eigenvalues`) + - Critical determinant |D_c| vs ψ_N (`plot_fixed_boundary_stability_criterion`) + - Least stable mode displacement (`plot_mode_displacement`) - Tearing stability Δ' at each rational surface (`plot_delta_prime`) -If no vacuum data is present (`vac_flag = false`), only the stability criterion and Δ' -panels are shown. +Panels the source cannot draw are left out; it throws an `ArgumentError` only when none can be drawn. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file + - `src`: A `ForceFreeStatesResult`, or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -552,28 +575,35 @@ panels are shown. A `Plots.jl` plot object. """ -function plot_ffs_summary(h5path; save_path=nothing) - has_vac = h5open(h5path, "r") do fid - haskey(fid, "ForceFreeStates/FreeBoundaryStability/W_freeboundary_eigenmodes") && - !isempty(read(fid["ForceFreeStates/FreeBoundaryStability/W_freeboundary_eigenmodes"])) - end - - p_crit = plot_fixed_boundary_stability_criterion(h5path) - p_dp = plot_delta_prime(h5path) - - if has_vac - p_evec = plot_energy_eigenvectors(h5path; matrix_type=:total) - p_modes = plot_mode_displacement(h5path) - p = plot(p_evec, p_crit, p_modes, p_dp; layout=(2, 2), size=(1100, 900), right_margin=10Plots.mm) - else - title!(p_crit, "Stability criterion (no vacuum data — rerun with vac_flag = true)") - p = plot(p_crit, p_dp; layout=(1, 2), size=(1100, 500), right_margin=10Plots.mm) - end - +function plot_ffs_summary(src::Source; save_path=nothing) + p = _summary(_ffs_panels(src), capabilities(src), plot_ffs_summary) isnothing(save_path) || savefig(p, save_path) return p end +# Draw the panels `caps` allows, two per row; throw naming every reason when none can be drawn. +function _summary(panels, caps, fn) + drawn = [draw() for (panel_fn, draw) in panels if _can(panel_fn, caps)] + if isempty(drawn) + missing = unique(reduce(vcat, (_missing(_entry(panel_fn), caps) for (panel_fn, _) in panels))) + throw(ArgumentError("$(nameof(fn)): no panel can be drawn: $(_why(missing))")) + end + n = length(drawn) + n <= 2 && return plot(drawn...; layout=(1, n), size=(550n, 500), right_margin=10Plots.mm) + n == 3 && push!(drawn, plot(; framestyle=:none)) + return plot(drawn...; layout=(2, 2), size=(1100, 900), right_margin=10Plots.mm) +end + +register!(plot_ffs_summary, (ForceFreeStatesResult, AbstractString), (), "the stability panels the run supplies") +register!(plot_mode_displacement, (ForceFreeStatesResult, AbstractString), (:solution_profiles,), "|ξ_ψ| per m of the least stable mode") +register!(plot_fixed_boundary_stability_criterion, (ForceFreeStatesResult, AbstractString), (:fixed_boundary_criterion,), "critical determinant against ψ_N") +register!(plot_energy_eigenvectors, (ForceFreeStatesResult, AbstractString), (:free_boundary,), "free-boundary energy eigenvectors per m") +register!(plot_eigenvalues, (ForceFreeStatesResult, AbstractString), (:free_boundary,), "free-boundary energy eigenvalues") +register!(plot_edge_stability_scan, (ForceFreeStatesResult, AbstractString), (:edge_scan,), "edge energies against the truncation q") +register!(plot_delta_prime, (ForceFreeStatesResult, AbstractString), (:delta_prime,), "Δ′_ii at each rational surface") +register!(plot_cond_fbar, (ForceFreeStatesResult, AbstractString), (:kinetic_scan,), "cond(F̄) scan for kinetic singular surfaces") +register!(plot_ballooning_alpha_boundary, (LocalStabilityResult, AbstractString), (:ballooning_scan,), "ballooning α against its first stability boundary") + """ plot_ballooning_alpha_boundaries(bnd; save_path=nothing, psi_min=0.0) diff --git a/src/Analysis/PerturbedEquilibrium.jl b/src/Analysis/PerturbedEquilibrium.jl index 83d5cde7b..76fd0aa85 100644 --- a/src/Analysis/PerturbedEquilibrium.jl +++ b/src/Analysis/PerturbedEquilibrium.jl @@ -1,8 +1,9 @@ """ PerturbedEquilibrium -Post-processing and visualization functions for GPEC perturbed equilibrium results stored -in the `PerturbedEquilibrium/` group of a GPEC HDF5 output file. +Post-processing and visualization functions for GPEC perturbed equilibrium results. Each plot takes the +`PerturbedEquilibriumState` in memory or the path of a `gpec.h5` it was saved to, and throws an `ArgumentError` +naming what to change in the run when the source lacks the data it draws (see `Analysis.available_plots`). """ module PerturbedEquilibrium @@ -10,26 +11,40 @@ using HDF5 using LaTeXStrings using Plots -# Check that a PE dataset exists and is non-empty. -function _has_pe_data(h5path, key) +using ..Analysis: register!, _require, _can, capabilities, REASONS +import ...PerturbedEquilibrium: PerturbedEquilibriumState + +const Source = Union{AbstractString,PerturbedEquilibriumState} + +const _COUPLING = "PerturbedEquilibrium/SingularCoupling/" +const _RESPONSE = "PerturbedEquilibrium/Response/" + +function _coupling_data(s::PerturbedEquilibriumState) + r = s.response + return (; field=s.resonant_area_weighted_field, current=s.resonant_current, island_hw=s.island_half_width, chirikov=s.chirikov_parameter, + delta_prime=s.delta_prime, psi=r.rational_psi, q=r.rational_q, n=r.rational_n) +end + +function _coupling_data(h5path::AbstractString) h5open(h5path, "r") do fid - haskey(fid, key) && !isempty(read(fid[key])) + rd(name) = read(fid[_COUPLING*name]) + (; field=rd("resonant_area_weighted_field"), current=rd("resonant_current"), island_hw=rd("island_half_width"), chirikov=rd("chirikov_parameter"), + delta_prime=rd("Delta_prime"), psi=rd("rational_psi"), q=rd("rational_q"), n=rd("rational_n")) end end """ - plot_resonant_area_weighted_field_amplitude(h5path; save_path=nothing) + plot_resonant_area_weighted_field_amplitude(src; save_path=nothing) Scatter plot of the resonant area-weighted field `|b^r|` per rational surface vs ψ_N. One marker series per toroidal mode n. Integer-valued q rational surfaces are annotated. The resonant area-weighted field `b^r = Φ^r/A^r` is the resonant flux normalized by the scalar surface area, in tesla [Pharr 2026]. -Requires the perturbed equilibrium module to have been run and -`SingularCoupling/resonant_area_weighted_field` to be present in the HDF5 file. +Requires the perturbed equilibrium singular coupling. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file with perturbed equilibrium output + - `src`: A `PerturbedEquilibriumState`, or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -39,17 +54,15 @@ Requires the perturbed equilibrium module to have been run and A `Plots.jl` plot object. """ -function plot_resonant_area_weighted_field_amplitude(h5path; save_path=nothing) - base = "PerturbedEquilibrium/SingularCoupling/" - _has_pe_data(h5path, base * "resonant_area_weighted_field") || - return plot(; title="No resonant area-weighted field data — run with perturbed equilibrium enabled", legend=false) - - resonant_area_weighted_field, rational_psi, rational_q, rational_n = h5open(h5path, "r") do fid - read(fid[base*"resonant_area_weighted_field"]), - read(fid[base*"rational_psi"]), - read(fid[base*"rational_q"]), - read(fid[base*"rational_n"]) - end +function plot_resonant_area_weighted_field_amplitude(src::Source; save_path=nothing) + _require(plot_resonant_area_weighted_field_amplitude, src) + p = _plot_resonant_area_weighted_field_amplitude(_coupling_data(src)) + isnothing(save_path) || savefig(p, save_path) + return p +end + +function _plot_resonant_area_weighted_field_amplitude(d) + resonant_area_weighted_field, rational_psi, rational_q, rational_n = d.field, d.psi, d.q, d.n p = plot(; xlabel="Norm. Poloidal Flux", ylabel="|b^r| [T]", title="Resonant area-weighted field |b^r| per surface", legend=:outertopright, @@ -65,22 +78,20 @@ function plot_resonant_area_weighted_field_amplitude(h5path; save_path=nothing) text(" q=$(round(Int, rational_q[k]))", 8, :left, :black)) end end - - isnothing(save_path) || savefig(p, save_path) return p end """ - plot_island_widths(h5path; save_path=nothing) + plot_island_widths(src; save_path=nothing) Scatter plot of island half-width `w/2` per rational surface vs ψ_N. Integer-valued q rational surfaces are annotated. -Requires `SingularCoupling/island_half_width` in the HDF5 file. +Requires the perturbed equilibrium singular coupling. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file with perturbed equilibrium output + - `src`: A `PerturbedEquilibriumState`, or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -90,16 +101,15 @@ Requires `SingularCoupling/island_half_width` in the HDF5 file. A `Plots.jl` plot object. """ -function plot_island_widths(h5path; save_path=nothing) - base = "PerturbedEquilibrium/SingularCoupling/" - _has_pe_data(h5path, base * "island_half_width") || - return plot(; title="No island width data — run with perturbed equilibrium enabled", legend=false) - - island_hw, rational_psi, rational_q = h5open(h5path, "r") do fid - read(fid[base*"island_half_width"]), - read(fid[base*"rational_psi"]), - read(fid[base*"rational_q"]) - end +function plot_island_widths(src::Source; save_path=nothing) + _require(plot_island_widths, src) + p = _plot_island_widths(_coupling_data(src)) + isnothing(save_path) || savefig(p, save_path) + return p +end + +function _plot_island_widths(d) + island_hw, rational_psi, rational_q = d.island_hw, d.psi, d.q p = scatter( rational_psi, island_hw; @@ -118,23 +128,21 @@ function plot_island_widths(h5path; save_path=nothing) annotate!(p, rational_psi[k], island_hw[k], text(" q=$(round(Int, rational_q[k]))", 8, :left, :black)) end - - isnothing(save_path) || savefig(p, save_path) return p end """ - plot_chirikov_parameter(h5path; save_path=nothing) + plot_chirikov_parameter(src; save_path=nothing) Scatter plot of the Chirikov overlap parameter per rational surface vs ψ_N, with a horizontal reference line at K = 1 (island overlap threshold). Points are colored red when K > 1. Integer-valued q rational surfaces are annotated. -Requires `SingularCoupling/chirikov_parameter` in the HDF5 file. +Requires the perturbed equilibrium singular coupling. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file with perturbed equilibrium output + - `src`: A `PerturbedEquilibriumState`, or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -144,16 +152,15 @@ Requires `SingularCoupling/chirikov_parameter` in the HDF5 file. A `Plots.jl` plot object. """ -function plot_chirikov_parameter(h5path; save_path=nothing) - base = "PerturbedEquilibrium/SingularCoupling/" - _has_pe_data(h5path, base * "chirikov_parameter") || - return plot(; title="No Chirikov data — run with perturbed equilibrium enabled", legend=false) - - chirikov, rational_psi, rational_q = h5open(h5path, "r") do fid - read(fid[base*"chirikov_parameter"]), - read(fid[base*"rational_psi"]), - read(fid[base*"rational_q"]) - end +function plot_chirikov_parameter(src::Source; save_path=nothing) + _require(plot_chirikov_parameter, src) + p = _plot_chirikov_parameter(_coupling_data(src)) + isnothing(save_path) || savefig(p, save_path) + return p +end + +function _plot_chirikov_parameter(d) + chirikov, rational_psi, rational_q = d.chirikov, d.psi, d.q colors = [k > 1.0 ? :red : :steelblue for k in chirikov] @@ -175,28 +182,24 @@ function plot_chirikov_parameter(h5path; save_path=nothing) annotate!(p, rational_psi[k], chirikov[k], text(" q=$(round(Int, rational_q[k]))", 8, :left, :black)) end - - isnothing(save_path) || savefig(p, save_path) return p end """ - plot_driven_delta_prime(h5path; save_path=nothing) + plot_driven_delta_prime(src; save_path=nothing) Scatter plot of `Re(Δ')` per rational surface vs ψ_N, computed by the perturbed -equilibrium module (from `PerturbedEquilibrium/SingularCoupling/Delta_prime`). +equilibrium singular coupling. One marker series per toroidal mode n. Integer-valued q rational surfaces are annotated. This is the forcing-driven Δ' (response to the applied perturbation amplitudes -in `intr.forcing_modes`); for the equilibrium-intrinsic Δ' from the STRIDE BVP, -read `SingularSurfaces/Delta_prime_matrix` from the HDF5 directly. - -Requires `PerturbedEquilibrium/SingularCoupling/Delta_prime` in the HDF5 file. +in `forcing_modes`); for the equilibrium-intrinsic Δ' of a Riccati or Galerkin solve, +see `ForceFreeStates.plot_delta_prime`. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file + - `src`: A `PerturbedEquilibriumState`, or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -206,17 +209,15 @@ Requires `PerturbedEquilibrium/SingularCoupling/Delta_prime` in the HDF5 file. A `Plots.jl` plot object. """ -function plot_driven_delta_prime(h5path; save_path=nothing) - base = "PerturbedEquilibrium/SingularCoupling/" - _has_pe_data(h5path, base * "Delta_prime") || - return plot(; title="No PE Δ' data — run with perturbed equilibrium enabled", legend=false) - - delta_prime, rational_psi, rational_q, rational_n = h5open(h5path, "r") do fid - read(fid[base*"Delta_prime"]), - read(fid[base*"rational_psi"]), - read(fid[base*"rational_q"]), - read(fid[base*"rational_n"]) - end +function plot_driven_delta_prime(src::Source; save_path=nothing) + _require(plot_driven_delta_prime, src) + p = _plot_driven_delta_prime(_coupling_data(src)) + isnothing(save_path) || savefig(p, save_path) + return p +end + +function _plot_driven_delta_prime(d) + delta_prime, rational_psi, rational_q, rational_n = d.delta_prime, d.psi, d.q, d.n p = plot(; xlabel="Norm. Poloidal Flux", ylabel="Re(Δ')", title="Tearing stability Δ' (driven, perturbed equilibrium)", legend=:outertopright, @@ -235,13 +236,11 @@ function plot_driven_delta_prime(h5path; save_path=nothing) text(" q=$(round(Int, rational_q[idx]))", 8, :left, :black)) end end - - isnothing(save_path) || savefig(p, save_path) return p end """ - plot_resonant_area_weighted_field(h5path; save_path=nothing) + plot_resonant_area_weighted_field(src; save_path=nothing) Five-panel summary of resonant coupling quantities at each rational surface vs ψ_N: @@ -255,7 +254,7 @@ Inspired by `plot_resonant_flux.py` from OMFIT GPEC. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file with perturbed equilibrium output + - `src`: A `PerturbedEquilibriumState`, or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -265,12 +264,14 @@ Inspired by `plot_resonant_flux.py` from OMFIT GPEC. A `Plots.jl` plot object. """ -function plot_resonant_area_weighted_field(h5path; save_path=nothing) - p1 = plot_resonant_area_weighted_field_amplitude(h5path) - p2 = plot_driven_delta_prime(h5path) - p3 = _plot_resonant_current(h5path) - p4 = plot_island_widths(h5path) - p5 = plot_chirikov_parameter(h5path) +function plot_resonant_area_weighted_field(src::Source; save_path=nothing) + _require(plot_resonant_area_weighted_field, src) + d = _coupling_data(src) + p1 = _plot_resonant_area_weighted_field_amplitude(d) + p2 = _plot_driven_delta_prime(d) + p3 = _plot_resonant_current(d) + p4 = _plot_island_widths(d) + p5 = _plot_chirikov_parameter(d) # Move outer legends inside so all panels have the same right margin for p in (p1, p2, p3) @@ -292,17 +293,8 @@ function plot_resonant_area_weighted_field(h5path; save_path=nothing) end # Internal helper — resonant current scatter plot -function _plot_resonant_current(h5path) - base = "PerturbedEquilibrium/SingularCoupling/" - _has_pe_data(h5path, base * "resonant_current") || - return plot(; title="No resonant current data", legend=false) - - resonant_current, rational_psi, rational_q, rational_n = h5open(h5path, "r") do fid - read(fid[base*"resonant_current"]), - read(fid[base*"rational_psi"]), - read(fid[base*"rational_q"]), - read(fid[base*"rational_n"]) - end +function _plot_resonant_current(d) + resonant_current, rational_psi, rational_q, rational_n = d.current, d.psi, d.q, d.n p = plot(; xlabel="Norm. Poloidal Flux", ylabel="|I_res|", title="Resonant current |I_res| per surface", legend=:outertopright, @@ -323,8 +315,39 @@ function _plot_resonant_current(h5path) return p end +const _COMPONENTS = Dict( + :xi_psi => (dataset="xi_psi", field=:psi, modes=:xi_modes, needs=:response_xi), + :b_psi => (dataset="b_psi_area_weighted", field=:b_psi_area_weighted, modes=:b_modes, needs=:response_b), # area-normalized b^ψ + :Jb_theta => (dataset="Jb_theta", field=:theta, modes=:b_modes, needs=:response_b), + :Jb_zeta => (dataset="Jb_zeta", field=:zeta, modes=:b_modes, needs=:response_b) +) + +function _response_data(s::PerturbedEquilibriumState, component) + c = _COMPONENTS[component] + ffs = s.response.ffs + return (; data=getfield(getfield(s, c.modes), c.field), psi=s.psi_grid, mlow=ffs.mlow, mhigh=ffs.mhigh, nhigh=ffs.nhigh, + q95=ffs.equil.params.q95, rational_psi=s.response.rational_psi) +end + +function _response_data(h5path::AbstractString, component) + h5open(h5path, "r") do fid + (; data=read(fid[_RESPONSE*_COMPONENTS[component].dataset]), psi=read(fid[_RESPONSE*"psi"]), mlow=read(fid["Info/mlow"]), + mhigh=read(fid["Info/mhigh"]), nhigh=read(fid["Info/nhigh"]), q95=read(fid["Equilibrium/q_95"]), + rational_psi=haskey(fid, _COUPLING * "rational_psi") ? read(fid[_COUPLING*"rational_psi"]) : Float64[]) + end +end + +# Check `component` names a response field `src` holds; `_require` covers ξ, which every component's run also has. +function _require_component(fn, src, component) + haskey(_COMPONENTS, component) || throw(ArgumentError("component must be one of :xi_psi, :b_psi, :Jb_theta, :Jb_zeta")) + _require(fn, src) + need = _COMPONENTS[component].needs + need in capabilities(src) || throw(ArgumentError("$(nameof(fn)): $(REASONS[need])")) + return nothing +end + """ - plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) + plot_mode_spectrogram(src; component=:xi_psi, save_path=nothing) Two-panel spectrogram of a perturbed equilibrium response field component: @@ -337,13 +360,13 @@ Inspired by `plot_spectrograms.py` from OMFIT GPEC. ### Arguments - - `h5path`: Path to a GPEC HDF5 output file with perturbed equilibrium response output + - `src`: A `PerturbedEquilibriumState` with the reconstructed response, or the path of a `gpec.h5` it was saved to ### Keyword arguments - `component`: Response field component to plot; one of `:xi_psi`, `:b_psi`, - `:Jb_theta`, `:Jb_zeta` (default: `:xi_psi`). `:b_psi` reads the area-normalized - `b_psi_area_weighted` dataset; `:Jb_theta` and `:Jb_zeta` are the Jacobian-weighted + `:Jb_theta`, `:Jb_zeta` (default: `:xi_psi`). `:b_psi` is the area-normalized + `b_psi_area_weighted` field; `:Jb_theta` and `:Jb_zeta` are the Jacobian-weighted contravariant field components. - `save_path`: If provided, save the figure to this path (default: `nothing`) @@ -351,36 +374,20 @@ Inspired by `plot_spectrograms.py` from OMFIT GPEC. A `Plots.jl` plot object. """ -function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) - comp_map = Dict( - :xi_psi => "xi_psi", - :b_psi => "b_psi_area_weighted", # area-normalized b^ψ - :Jb_theta => "Jb_theta", - :Jb_zeta => "Jb_zeta" - ) - haskey(comp_map, component) || - error("component must be one of :xi_psi, :b_psi, :Jb_theta, :Jb_zeta") - - base = "PerturbedEquilibrium/Response/" - dataset_path = base * comp_map[component] - - _has_pe_data(h5path, dataset_path) || - return plot(; title="No response data — run with perturbed equilibrium enabled", legend=false) - - data, psi_response, mlow, mhigh, nhigh, q95, rational_psi = h5open(h5path, "r") do fid - read(fid[dataset_path]), - read(fid[base*"psi"]), - read(fid["Info/mlow"]), read(fid["Info/mhigh"]), read(fid["Info/nhigh"]), - read(fid["Equilibrium/q_95"]), - read(fid["PerturbedEquilibrium/SingularCoupling/rational_psi"]) - end +function plot_mode_spectrogram(src::Source; component=:xi_psi, save_path=nothing) + _require_component(plot_mode_spectrogram, src, component) + p = _plot_mode_spectrogram(_response_data(src, component), component) + isnothing(save_path) || savefig(p, save_path) + return p +end - mpert = mhigh - mlow + 1 - m_vals = mlow:mhigh - data_mn = data[:, 1:mpert] # (npsi, mpert): first mpert columns = first n's modes +function _plot_mode_spectrogram(d, component) + mpert = d.mhigh - d.mlow + 1 + m_vals = d.mlow:d.mhigh + data_mn = d.data[:, 1:mpert] # (npsi, mpert): first mpert columns = first n's modes # Top panel: line plot per mode — only label resonant range m ∈ [0, nhigh·q95) - m_max_legend = nhigh * q95 + m_max_legend = d.nhigh * d.q95 p1 = plot(; xlabel="Norm. Poloidal Flux", ylabel="|$(component)|", @@ -392,13 +399,13 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) cmap = cgrad(:roma, mpert; categorical=true) for (i, m) in enumerate(m_vals) show_label = 0 <= m < m_max_legend - plot!(p1, psi_response, abs.(data_mn[:, i]); + plot!(p1, d.psi, abs.(data_mn[:, i]); label=(show_label ? "m=$m" : nothing), color=cmap[i], linewidth=1.5) end # Bottom panel: heatmap in (m, ψ_N) space — psi on vertical axis, m on horizontal p2 = heatmap( - collect(m_vals), psi_response, abs.(data_mn); + collect(m_vals), d.psi, abs.(data_mn); xlabel="m", ylabel="Norm. Poloidal Flux", title="", @@ -408,28 +415,27 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) bottom_margin=5Plots.mm ) # Overlay rational surface locations as white dashed lines - for psi_s in rational_psi + for psi_s in d.rational_psi hline!(p2, [psi_s]; linestyle=:dash, color=:white, linewidth=1.5, label=nothing) end - p = plot(p1, p2; layout=(2, 1), size=(950, 750)) - - isnothing(save_path) || savefig(p, save_path) - return p + return plot(p1, p2; layout=(2, 1), size=(950, 750)) end """ - plot_perturbed_equilibrium_summary(h5path; save_path=nothing) + plot_perturbed_equilibrium_summary(src; save_path=nothing) -Three-panel composite summary of perturbed equilibrium results: +Composite summary of perturbed equilibrium results, with each panel the source can supply: - Top-left: Resonant area-weighted field per surface (`plot_resonant_area_weighted_field_amplitude`) - Top-right: Edge |b_ψ| spectrum - Bottom: ξ_ψ mode spectrogram (`plot_mode_spectrogram`) +Panels the source cannot draw are left out; it throws an `ArgumentError` only when none can be drawn. + ### Arguments - - `h5path`: Path to a GPEC HDF5 output file with perturbed equilibrium output + - `src`: A `PerturbedEquilibriumState`, or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -439,32 +445,31 @@ Three-panel composite summary of perturbed equilibrium results: A `Plots.jl` plot object. """ -function plot_perturbed_equilibrium_summary(h5path; save_path=nothing) - p_islands = plot_resonant_area_weighted_field_amplitude(h5path) - p_bpsi = _plot_bpsi_edge_spectrum(h5path) - p_spectro = plot_mode_spectrogram(h5path; component=:xi_psi) - - l = @layout [grid(1, 2){0.35h}; b] - p = plot(p_islands, p_bpsi, p_spectro; layout=l, size=(1100, 1100)) +function plot_perturbed_equilibrium_summary(src::Source; save_path=nothing) + caps = capabilities(src) + top = Plots.Plot[] + _can(plot_resonant_area_weighted_field_amplitude, caps) && push!(top, _plot_resonant_area_weighted_field_amplitude(_coupling_data(src))) + :response_b in caps && push!(top, _plot_bpsi_edge_spectrum(_response_data(src, :b_psi))) + has_spectro = _can(plot_mode_spectrogram, caps) + if !has_spectro && isempty(top) + throw(ArgumentError("plot_perturbed_equilibrium_summary: no panel can be drawn: $(REASONS[:singular_coupling]); $(REASONS[:response_xi])")) + end + p = if !has_spectro + plot(top...; layout=(1, length(top)), size=(550 * length(top), 450)) + else + p_spectro = _plot_mode_spectrogram(_response_data(src, :xi_psi), :xi_psi) + isempty(top) ? p_spectro : plot(top..., p_spectro; layout=@layout([grid(1, length(top)){0.35h}; b]), size=(1100, 1100)) + end isnothing(save_path) || savefig(p, save_path) return p end # Internal helper — |b_psi(m)| spectrum at the outermost psi surface -function _plot_bpsi_edge_spectrum(h5path) - base = "PerturbedEquilibrium/Response/" - _has_pe_data(h5path, base * "b_psi_area_weighted") || - return plot(; title="No b_psi data — run with perturbed equilibrium enabled", legend=false) - - data, mlow, mhigh = h5open(h5path, "r") do fid - read(fid[base*"b_psi_area_weighted"]), - read(fid["Info/mlow"]), read(fid["Info/mhigh"]) - end - - mpert = mhigh - mlow + 1 - m_vals = mlow:mhigh - bpsi_edge = abs.(data[end, 1:mpert]) +function _plot_bpsi_edge_spectrum(d) + mpert = d.mhigh - d.mlow + 1 + m_vals = d.mlow:d.mhigh + bpsi_edge = abs.(d.data[end, 1:mpert]) p = bar( collect(m_vals), bpsi_edge; @@ -481,4 +486,13 @@ function _plot_bpsi_edge_spectrum(h5path) return p end +const _PE_SOURCES = (PerturbedEquilibriumState, AbstractString) +register!(plot_perturbed_equilibrium_summary, _PE_SOURCES, (), "resonant field, edge b_ψ spectrum and ξ spectrogram") +register!(plot_resonant_area_weighted_field, _PE_SOURCES, (:singular_coupling,), "five resonant-coupling panels per rational surface") +register!(plot_resonant_area_weighted_field_amplitude, _PE_SOURCES, (:singular_coupling,), "|b̄_res| per rational surface") +register!(plot_island_widths, _PE_SOURCES, (:singular_coupling,), "island half-width per rational surface") +register!(plot_chirikov_parameter, _PE_SOURCES, (:singular_coupling,), "Chirikov parameter per rational surface") +register!(plot_driven_delta_prime, _PE_SOURCES, (:singular_coupling,), "forcing-driven Δ′ per rational surface") +register!(plot_mode_spectrogram, _PE_SOURCES, (:response_xi,), "|response| per m against ψ_N (component = :xi_psi, :b_psi, :Jb_theta, :Jb_zeta)") + end # module PerturbedEquilibrium diff --git a/src/Analysis/Registry.jl b/src/Analysis/Registry.jl new file mode 100644 index 000000000..4f72964ac --- /dev/null +++ b/src/Analysis/Registry.jl @@ -0,0 +1,231 @@ +# Registry of the Analysis plots: which data products each one draws, which sources (a result in memory or a +# gpec.h5 path) supply them, and what a source must have computed for the plot to draw. Each plotting submodule +# registers its own plots at load time; `capabilities` is the single record of what each source holds. + +using HDF5 + +import ..Equilibrium: PlasmaEquilibrium +import ..LocalStability: LocalStabilityResult +import ..ForceFreeStates: ForceFreeStatesResult +import ..PerturbedEquilibrium: PerturbedEquilibriumState +import ..ErrorFields as EF + +""" + PlotEntry + +A registered plot. + +## Fields + + - `fn::Function` - The plotting function. + - `sources::Tuple` - The types it accepts as its first argument; `AbstractString` means a `gpec.h5` path. + - `needs::Tuple{Vararg{Symbol}}` - The capabilities (see [`capabilities`](@ref)) a source must have for it to draw. + - `description::String` - One line on what it draws. +""" +struct PlotEntry + fn::Function + sources::Tuple + needs::Tuple{Vararg{Symbol}} + description::String +end + +const PLOTS = PlotEntry[] + +# Add a plot to the registry; called by each plotting submodule as it loads. +register!(fn, sources, needs, description) = (push!(PLOTS, PlotEntry(fn, Tuple(sources), Tuple(needs), description)); fn) + +""" +What each capability needs, worded as what to change in the run to get it. +""" +const REASONS = Dict{Symbol,String}( + :equilibrium => "needs an equilibrium", + :singular_surfaces => "needs rational surfaces in the integration domain", + :solution_profiles => "needs stored ξ profiles: use the Forward integrator or a matched Galerkin solve (Riccati and unmatched Galerkin store none)", + :fixed_boundary_criterion => "needs the critical determinant scan: use the Forward or Riccati integrator (Galerkin computes none)", + :free_boundary => "needs the free-boundary energies: run with vac_flag = true and the Forward or Riccati integrator", + :edge_scan => "needs the edge stability scan: set dW_edge_scan_start below psilim", + :delta_prime => "needs the Δ′ matrix: use the Riccati or Galerkin integrator (Forward computes none)", + :kinetic_scan => "needs the kinetic singular-surface scan: run with a kinetic source and singfac_min > 0", + :ballooning_scan => "needs the ballooning scan: add a [LocalStability] section", + :gse_files => "needs gse.h5 next to the file: set diagnose_src = true in [Equilibrium]", + :singular_coupling => "needs the singular coupling: run PerturbedEquilibrium with compute_singular_coupling = true", + :response_xi => "needs the reconstructed displacement ξ: run PerturbedEquilibrium on a Forward free-boundary result", + :response_b => "needs the reconstructed perturbed field b: run PerturbedEquilibrium on a Forward free-boundary result", + :coil_sets => "needs coil forcing: give ForcingTerms coil sets (forcing_data_format = \"coil\")", + :forcing_modes => "needs forcing modes: give ForcingTerms a forcing", + :ef_sensitivities => "needs the coil sensitivities: add an [ErrorFields] section", + :dominant_mode => "needs the dominant resonant-coupling mode: run PerturbedEquilibrium with compute_singular_coupling = true", + :ef_monte_carlo => "needs the tolerance Monte Carlo: give ErrorFields a tolerance file", + :ef_risk => "needs the locking risk: give ErrorFields a tolerance file and a scenario", + :ef_scan => "needs the tolerance scan: set RiskControl scan_scales" +) + +""" + capabilities(src) -> Set{Symbol} + +The data products `src` holds, by the names plots declare in their `needs`: a result in memory, or a +`gpec.h5` path (a dataset counts when present and non-empty). See [`REASONS`](@ref) for what each name means. +""" +function capabilities end + +# A type no plot takes holds nothing a plot reads. +capabilities(_) = Set{Symbol}() + +capabilities(::PlasmaEquilibrium) = Set([:equilibrium]) + +function capabilities(r::ForceFreeStatesResult) + caps = Set([:equilibrium]) + isempty(r.surfaces) || push!(caps, :singular_surfaces) + r.solution === nothing || push!(caps, :solution_profiles) + diag = r.diagnostics + diag === nothing || isempty(diag.crit_store) || push!(caps, :fixed_boundary_criterion) + diag === nothing || isempty(diag.edge_scan.psi) || push!(caps, :edge_scan) + r.free_boundary === nothing || isempty(r.free_boundary.et) || push!(caps, :free_boundary) + r.delta_prime === nothing || isempty(r.delta_prime.matrix) || push!(caps, :delta_prime) + isempty(r.kinetic.scan_psi) || push!(caps, :kinetic_scan) + return caps +end + +capabilities(r::LocalStabilityResult) = isempty(r.alpha) ? Set{Symbol}() : Set([:ballooning_scan]) + +function capabilities(s::PerturbedEquilibriumState) + caps = Set{Symbol}() + isempty(s.resonant_area_weighted_field) || push!(caps, :singular_coupling) + s.xi_modes === nothing || push!(caps, :response_xi) + s.b_modes === nothing || push!(caps, :response_b) + isempty(s.coil_sets) || push!(caps, :coil_sets) + isempty(s.forcing_modes) || push!(caps, :forcing_modes) + return caps +end + +function capabilities(r::EF.ErrorFieldsResult) + caps = Set([:ef_sensitivities]) + isempty(r.response.dominant_singular_values) || push!(caps, :dominant_mode) + r.monte_carlo === nothing || push!(caps, :ef_monte_carlo) + r.risk === nothing || push!(caps, :ef_risk) + r.scan === nothing || push!(caps, :ef_scan) + return caps +end + +const _H5_CAPABILITIES = ( + :equilibrium => "Equilibrium/Profiles/psi", + :singular_surfaces => "SingularSurfaces/rational_psi", + :solution_profiles => ("ForceFreeStates/Solutions/ForwardIntegration/xi_psi", "ForceFreeStates/Solutions/GalerkinIntegration/xi_psi"), + :fixed_boundary_criterion => "ForceFreeStates/Solutions/ForwardIntegration/crit", + :edge_scan => "ForceFreeStates/EdgeScan/psi", + :free_boundary => "ForceFreeStates/FreeBoundaryStability/eigenmode_energies", + :delta_prime => "SingularSurfaces/Delta_prime_matrix", + :kinetic_scan => "SingularSurfaces/Kinetic/scan_psi", + :ballooning_scan => "LocalStability/alpha", + :singular_coupling => "PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field", + :response_xi => "PerturbedEquilibrium/Response/xi_psi", + :response_b => "PerturbedEquilibrium/Response/b_psi_area_weighted", + :ef_sensitivities => EF._H5_GROUP, + :dominant_mode => "PerturbedEquilibrium/SingularCoupling/DominantMode", + :ef_monte_carlo => EF._MC_GROUP, + :ef_risk => EF._RISK_GROUP, + :ef_scan => "$(EF._RISK_GROUP)/ToleranceScan" +) + +function capabilities(h5path::AbstractString) + isfile(h5path) || throw(ArgumentError("no such file: $h5path")) + caps = Set{Symbol}() + h5open(h5path, "r") do f + # A group counts when present; a dataset when it also holds data. + has(k) = haskey(f, k) && (f[k] isa HDF5.Group || !isempty(f[k])) + for (cap, keys) in _H5_CAPABILITIES + any(has, keys isa Tuple ? keys : (keys,)) && push!(caps, cap) + end + end + isfile(joinpath(dirname(h5path), "gse.h5")) && push!(caps, :gse_files) + return caps +end + +# The registered entry of `fn`. +function _entry(fn) + i = findfirst(e -> e.fn === fn, PLOTS) + i === nothing && error("$(nameof(fn)) is not a registered plot") + return PLOTS[i] +end + +_accepts(e::PlotEntry, src) = any(T -> src isa T, e.sources) +_missing(e::PlotEntry, caps) = [k for k in e.needs if k ∉ caps] +_why(missing) = join((REASONS[k] for k in missing), "; ") + +""" + _require(fn, src) + +Throw an `ArgumentError` naming what the run must change when `src` lacks a capability that plot `fn` needs. +""" +function _require(fn, src) + e = _entry(fn) + missing = _missing(e, capabilities(src)) + isempty(missing) || throw(ArgumentError("$(nameof(fn)): $(_why(missing))")) + return nothing +end + +# Whether `src` has every capability plot `fn` needs; composites use it to skip the panels they cannot draw. +_can(fn, caps) = isempty(_missing(_entry(fn), caps)) + +""" + PlotCatalog + +The plots registered for a source, from [`available_plots`](@ref): one section per result, each listing the plots +that accept it, paired with the capabilities it lacks for them (none when it can draw them). +""" +struct PlotCatalog + sections::Vector{Tuple{String,Vector{Tuple{PlotEntry,Vector{Symbol}}}}} +end + +_source_label(h5path::AbstractString) = h5path +_source_label(src) = string(nameof(typeof(src))) +_source_label(r::ForceFreeStatesResult) = "ForceFreeStatesResult ($(r.integrator))" + +function _section(src) + caps = capabilities(src) + rows = [(e, _missing(e, caps)) for e in PLOTS if _accepts(e, src)] + return (_source_label(src), rows) +end + +""" + available_plots(src) -> PlotCatalog + +The registered plots for `src` with whether each can draw from it, and the reason when it cannot. `src` is a +result in memory, a `gpec.h5` path, or the named tuple `main` returns (one section per stage result). Every +listed function takes `src` as its first argument. +""" +available_plots(src) = PlotCatalog([_section(src)]) +available_plots(results::NamedTuple) = + PlotCatalog([_section(x) for x in unique(values(results)) if x !== nothing && any(e -> _accepts(e, x), PLOTS)]) + +""" + plot_counts(src) -> (available, total) + +How many registered plots accept `src`, and how many of those it can draw. +""" +function plot_counts(src) + _, rows = _section(src) + return count(r -> isempty(r[2]), rows), length(rows) +end + +Base.length(c::PlotCatalog) = sum(length(rows) for (_, rows) in c.sections; init=0) + +_qualified(e::PlotEntry) = string(nameof(parentmodule(e.fn)), ".", nameof(e.fn)) + +function Base.show(io::IO, ::MIME"text/plain", c::PlotCatalog) + isempty(c.sections) && return print(io, "PlotCatalog: no registered plots accept this source") + for (i, (label, rows)) in enumerate(c.sections) + i > 1 && print(io, "\n") + ok = count(r -> isempty(r[2]), rows) + print(io, "Plots for ", label, ": ", ok, " of ", length(rows), " available") + width = maximum((length(_qualified(e)) for (e, _) in rows); init=0) + for (e, missing) in rows + mark = isempty(missing) ? "✓" : "✗" + note = isempty(missing) ? e.description : _why(missing) + print(io, "\n ", mark, " ", rpad(_qualified(e), width), " ", note) + end + end + return nothing +end + +Base.show(io::IO, c::PlotCatalog) = print(io, "PlotCatalog(", length(c), " plots)") diff --git a/src/Display.jl b/src/Display.jl index fff5181ad..28264f92a 100644 --- a/src/Display.jl +++ b/src/Display.jl @@ -1,7 +1,8 @@ # Concise REPL display for the published result types: a one-line `show(io, x)` used inside # containers (the named tuple `main` returns, nested fields) and a short labelled `text/plain` # summary. Arrays are summarized by shape or range and never printed. Absent products read -# "not computed" (not requested or not producible) or "none" (computed, but empty). +# "not computed" (not requested or not producible) or "none" (computed, but empty). Results with +# registered plots end with how many of them they can draw. # # Included directly into the GeneralizedPerturbedEquilibrium module, after every submodule. @@ -34,6 +35,13 @@ end _row(io, label, value) = print(io, "\n ", rpad(label, 19), value) +# How many registered plots `x` can draw, pointing at the catalog that names them. +function _plots_row(io, x) + ok, n = Analysis.plot_counts(x) + n == 0 || _row(io, "Plots", "$ok of $n available (see available_plots)") + return nothing +end + # Range of `ψ` where `mask` holds, as "k/N points, ψ ∈ [a, b]"; "none" when nowhere. function _where(psi, mask) k = count(mask) @@ -72,6 +80,7 @@ function Base.show(io::IO, ::MIME"text/plain", eq::PlasmaEquilibrium) _row(io, "Shape", "R = $(_fmt(p.rmean)) m, a = $(_fmt(p.amean)) m, κ = $(_fmt(p.kappa)), δ_upper = $(_fmt(p.delta1)), δ_lower = $(_fmt(p.delta2))") _row(io, "Grid", "$(length(eq.rzphi_xs)) ψ × $(length(eq.rzphi_ys)) θ, ψ ∈ $(_span(eq.rzphi_xs)), $(eq.config.jac_type) coordinates") _row(io, "Kinetic", eq.kinetic === nothing ? "none attached" : _list(_kinetic_names(eq.kinetic); max=14)) + _plots_row(io, eq) return nothing end @@ -103,6 +112,7 @@ function Base.show(io::IO, ::MIME"text/plain", r::LocalStability.LocalStabilityR _row(io, "Mercier", "D_I > 0 (unstable): " * _where(r.psi, r.D_I .> 0)) _row(io, "Resistive", "D_R > 0 (unstable): " * _where(r.psi, r.D_R .> 0)) _row(io, "Ballooning", _ballooning_note(r)) + _plots_row(io, r) return nothing end @@ -218,6 +228,7 @@ function Base.show(io::IO, ::MIME"text/plain", r::ForceFreeStates.ForceFreeState r.match === nothing || _row(io, "Match", sprint(show, r.match)) r.kinetic.kmsing > 0 && _row(io, "Kinetic", "$(_plural(r.kinetic.kmsing, "kinetic singular surface"))") _row(io, "Wall", _wall_label(r.wall_settings)) + _plots_row(io, r) return nothing end @@ -302,6 +313,7 @@ function Base.show(io::IO, ::MIME"text/plain", s::PerturbedEquilibriumState) _row(io, "Profiles", s.xi_modes === nothing ? "not computed" : "ξ, b on $(length(s.psi_grid)) ψ nodes") s.torque_response === nothing || _row(io, "Torque resp.", "computed") iszero(s.regularization_width) || _row(io, "Regularized", "width $(_fmt(s.regularization_width))") + _plots_row(io, s) return nothing end @@ -381,6 +393,7 @@ function Base.show(io::IO, ::MIME"text/plain", r::ErrorFieldsResult) end r.scan === nothing || _row(io, "Tolerance scan", "$(length(r.scan.tolerance_scale)) points, $(_in("scale", r.scan.tolerance_scale))") r.couplings === nothing || _row(io, "EFC coils", _list([c.coil_name for c in r.couplings])) + _plots_row(io, r) return nothing end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index bda2b7c8d..26db532d1 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -81,6 +81,7 @@ using .ForcingTerms: RMPField using .PerturbedEquilibrium: plasma_response, perturbed_equilibrium, PlasmaResponse using .KineticForces: kinetic_torque, KineticForcesControl, KineticForcesResult using .ErrorFields: error_fields, ErrorFieldsResult +using .Analysis: available_plots using .ForceFreeStates: ForceBalanceProblem, kinetic_matrices, _prepare_solve, maybe_reform_equilibrium, ForceFreeStatesControl, kinetic_source include("Deck.jl") @@ -308,6 +309,7 @@ export PlasmaEquilibrium, ForceBalanceProblem, Forward, Riccati, Galerkin, TestK export MatchProblem, TearingProblem, TearingConfig, TearingResult, GGJModel, SLAYERModel export LocalStabilityProblem, Ballooning, LocalStabilityResult export KineticForcesControl, KineticForcesResult, ErrorFieldsResult +export available_plots include("Precompile.jl") diff --git a/test/display_check.jl b/test/display_check.jl index 19cb89f3a..559bcf9ac 100644 --- a/test/display_check.jl +++ b/test/display_check.jl @@ -16,7 +16,9 @@ function _check_main_display(ret) for x in ret x === nothing && continue @test parentmodule(which(show, Tuple{IO,MIME"text/plain",typeof(x)})) === GPEC - _check_display(x) + plain = _check_display(x) + # A result with registered plots says how many it can draw. + @test occursin("Plots", plain) == (last(GPEC.Analysis.plot_counts(x)) > 0) end @test !occursin('\n', repr(ret)) end diff --git a/test/runtests_analysis_plots.jl b/test/runtests_analysis_plots.jl index 4739fe5ba..b2df0da29 100644 --- a/test/runtests_analysis_plots.jl +++ b/test/runtests_analysis_plots.jl @@ -2,15 +2,15 @@ using HDF5 using Plots using TOML -# Smoke tests of the gpec.h5 plots against real writer output: each plot must read the datasets the -# writers emit, so a renamed path fails here instead of silently drawing a "No ... data" placeholder. -@testset "Analysis plots from gpec.h5" begin +# Every registered plot against real writer output, from the result in memory and from the gpec.h5 it was saved to: +# both sources must supply the same data, draw the same plots, and refuse the same ones with the registry's reason. +@testset "Analysis plots from a result or gpec.h5" begin ENV["GKSwstype"] = "100" GPEC = GeneralizedPerturbedEquilibrium A = GPEC.Analysis - # Copy a coarse fixture deck into `dir` with HDF5 output on and the given extra sections; return its gpec.h5. - function _run_fixture(dir, fixture; sections=Dict{String,Any}(), ffs=Dict{String,Any}()) + # Copy a coarse fixture deck into `dir` with HDF5 output on and the given extra sections; return main's results and its gpec.h5. + function _run_fixture(dir, fixture; sections=Dict{String,Any}(), ffs=Dict{String,Any}(), edit=identity) template = joinpath(@__DIR__, "test_data", fixture) for name in readdir(template) cp(joinpath(template, name), joinpath(dir, name)) @@ -19,16 +19,58 @@ using TOML inputs = TOML.parsefile(toml_path) merge!(inputs, sections) merge!(inputs["ForceFreeStates"], ffs) + edit(inputs) open(io -> TOML.print(io, inputs), toml_path, "w") - GPEC.main([dir]) - return joinpath(dir, "gpec.h5") + return GPEC.main([dir]), joinpath(dir, "gpec.h5") end - _placeholder(p) = any(sp -> startswith(string(get(sp.attr, :title, "")), "No "), p.subplots) - _drew(p) = p isa Plots.Plot && !_placeholder(p) - @testset "ideal run with coil forcing and the perturbed-equilibrium response" begin + # The data extractors each dual-source plot reads, by the capability that guards them. + extractors = ( + (A.ForceFreeStates._mode_displacement_data, :solution_profiles), (A.ForceFreeStates._criterion_data, :fixed_boundary_criterion), + (A.ForceFreeStates._energy_eigenvector_data, :free_boundary), (A.ForceFreeStates._delta_prime_data, :delta_prime), + (A.ForceFreeStates._cond_fbar_data, :kinetic_scan), (A.Equilibrium._global_data, :equilibrium), + (A.PerturbedEquilibrium._coupling_data, :singular_coupling)) + _same(a::NamedTuple, b::NamedTuple) = keys(a) == keys(b) && all(_same(a[k], b[k]) for k in keys(a)) + _same(a::AbstractArray{<:Number}, b::AbstractArray{<:Number}) = size(a) == size(b) && isapprox(a, b; rtol=1e-10, nans=true) + _same(a::Number, b::Number) = isapprox(a, b; rtol=1e-12) + _same(a, b) = isequal(a, b) + + # Every plot registered for `x` has a method for each source it lists, draws from `x` and from `h5` when the + # capabilities allow, and otherwise throws the registry's reason from both. + function _check_sources(x, h5) + caps = A.capabilities(x) + # The file holds every capability a plot can read from a path (coil geometry is in memory only). + path_needs = union((e.needs for e in A.PLOTS if A._accepts(e, h5))...) + @test intersect(caps, path_needs) ⊆ A.capabilities(h5) + for (f, cap) in extractors + cap in caps && hasmethod(f, Tuple{typeof(x)}) && @test _same(f(x), f(h5)) + end + for e in A.PLOTS + A._accepts(e, x) || continue + @test all(T -> hasmethod(e.fn, Tuple{T}), e.sources) + sources = A._accepts(e, h5) ? (x, h5) : (x,) + if A._can(e.fn, caps) + for src in sources + @test e.fn(src) isa Plots.Plot + end + else + reason = A._why(A._missing(e, caps)) + for src in sources + err = try + e.fn(src) + nothing + catch caught + caught + end + @test err isa ArgumentError && occursin(reason, err.msg) + end + end + end + end + + @testset "ideal Forward run with coil forcing, the response and local stability" begin mktempdir() do dir - h5 = _run_fixture( + ret, h5 = _run_fixture( dir, "regression_tj_ideal_example"; ffs=Dict{String,Any}("dW_edge_scan_start" => 0.95), @@ -38,34 +80,31 @@ using TOML "currents" => [2.0e3], "tiltx" => [3.0])]), "PerturbedEquilibrium" => Dict{String,Any}("compute_singular_coupling" => true)) ) - - for f in (A.Equilibrium.plot_qprofile, A.Equilibrium.plot_pressure_profile, A.Equilibrium.plot_f_profile, - A.Equilibrium.plot_flux_surfaces, A.Equilibrium.plot_equilibrium_summary) - @test _drew(f(h5)) - end - for f in (A.ForceFreeStates.plot_mode_displacement, A.ForceFreeStates.plot_fixed_boundary_stability_criterion, - A.ForceFreeStates.plot_energy_eigenvectors, A.ForceFreeStates.plot_edge_stability_scan, - A.ForceFreeStates.plot_eigenvalues, A.ForceFreeStates.plot_delta_prime, - A.ForceFreeStates.plot_ballooning_alpha_boundary, A.ForceFreeStates.plot_ffs_summary) - @test _drew(f(h5)) - end - for f in (A.PerturbedEquilibrium.plot_resonant_area_weighted_field_amplitude, A.PerturbedEquilibrium.plot_island_widths, - A.PerturbedEquilibrium.plot_chirikov_parameter, A.PerturbedEquilibrium.plot_driven_delta_prime, - A.PerturbedEquilibrium.plot_resonant_area_weighted_field, A.PerturbedEquilibrium.plot_perturbed_equilibrium_summary) - @test _drew(f(h5)) + for x in (ret.equil, ret.ffs, ret.pe) + _check_sources(x, h5) end + ret.local_stability === nothing || _check_sources(ret.local_stability, h5) for component in (:xi_psi, :b_psi, :Jb_theta, :Jb_zeta) - @test _drew(A.PerturbedEquilibrium.plot_mode_spectrogram(h5; component=component)) + @test _same(A.PerturbedEquilibrium._response_data(ret.pe, component), A.PerturbedEquilibrium._response_data(h5, component)) + @test A.PerturbedEquilibrium.plot_mode_spectrogram(ret.pe; component) isa Plots.Plot end - @test_throws ErrorException A.PerturbedEquilibrium.plot_mode_spectrogram(h5; component=:b_theta) + @test_throws ArgumentError A.PerturbedEquilibrium.plot_mode_spectrogram(h5; component=:b_theta) - # Optional outputs this run does not produce: no kinetic scan, no GSE sidecar files. - @test (@test_logs (:warn,) A.ForceFreeStates.plot_cond_fbar(h5)) === nothing - @test A.Equilibrium.plot_gse_by_theta(h5) === nothing - @test A.Equilibrium.plot_gse_integrated(h5) === nothing + # The catalog lists a Forward run's absent products with what to change, for a result, a path and main's tuple. + listing = repr(MIME"text/plain"(), available_plots(ret.ffs)) + @test occursin("✓ ForceFreeStates.plot_mode_displacement", listing) + @test occursin("✗ ForceFreeStates.plot_delta_prime", listing) && occursin(A.REASONS[:delta_prime], listing) + @test occursin("✗ Equilibrium.plot_gse_by_theta", repr(MIME"text/plain"(), available_plots(h5))) + @test length(available_plots(ret).sections) == count(x -> x !== nothing && any(e -> A._accepts(e, x), A.PLOTS), unique(values(ret))) + @test A.plot_counts(ret.ffs) == (11, 13) + + # Optional outputs this run does not produce throw instead of drawing a placeholder. + @test_throws ArgumentError A.ForceFreeStates.plot_cond_fbar(h5) + @test_throws ArgumentError A.Equilibrium.plot_gse_by_theta(h5) + @test_throws ArgumentError A.Equilibrium.plot_gse_integrated(h5) png = joinpath(dir, "qprofile.png") - A.Equilibrium.plot_qprofile(h5; save_path=png) + A.Equilibrium.plot_qprofile(ret.equil; save_path=png) @test isfile(png) theta_data, theta_grid, n_vals = A.PerturbedEquilibriumModes.modes_to_theta(h5, "PerturbedEquilibrium/Response/xi_psi") @@ -74,6 +113,19 @@ using TOML end end + @testset "Riccati run offers Δ′ but no solution profiles" begin + mktempdir() do dir + ret, h5 = _run_fixture(dir, "regression_solovev_ideal_example"; + edit=d -> (d["ForceFreeStates"]["integrator"] = "Riccati"; d["Riccati"] = pop!(d, "Forward"))) + _check_sources(ret.ffs, h5) + @test :delta_prime in A.capabilities(ret.ffs) && :solution_profiles ∉ A.capabilities(ret.ffs) + @test_throws A.REASONS[:solution_profiles] A.ForceFreeStates.plot_mode_displacement(ret.ffs) + @test occursin("✗ ForceFreeStates.plot_mode_displacement", repr(MIME"text/plain"(), available_plots(ret.ffs))) + # The summary leaves out the ξ panel rather than failing: three panels and a blank fourth. + @test length(A.ForceFreeStates.plot_ffs_summary(ret.ffs).subplots) == 4 + end + end + @testset "kinetic run draws the F̄ condition-number scan" begin # Synthetic test kinetic matrices are script-only, so this run goes through the API and saves its result. mktempdir() do dir @@ -81,9 +133,11 @@ using TOML fixture = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") kwargs = merge(GPEC.force_balance_keywords(inputs, fixture), Dict(:kinetic => TestKineticMatrices(; scale=1e-9))) prob = ForceBalanceProblem(PlasmaEquilibrium(config, analytic); kwargs...) + ffs = solve(prob, GPEC.force_balance_algorithm(inputs)) h5 = joinpath(dir, "gpec.h5") - save_to_hdf5(h5, solve(prob, GPEC.force_balance_algorithm(inputs))) - @test _drew(A.ForceFreeStates.plot_cond_fbar(h5)) + save_to_hdf5(h5, ffs) + _check_sources(ffs, h5) + @test A.ForceFreeStates.plot_cond_fbar(ffs) isa Plots.Plot end end end diff --git a/test/runtests_error_fields.jl b/test/runtests_error_fields.jl index 35f01a3c7..af15828d0 100644 --- a/test/runtests_error_fields.jl +++ b/test/runtests_error_fields.jl @@ -259,6 +259,20 @@ include("display_check.jl") @test AEF.plot_coil_sensitivities(h5path; quantity=:rim)[1][:yaxis][:guide] == "|δ| per mm of rim displacement" @test_throws ArgumentError AEF.plot_coil_sensitivities(h5path; quantity=:bogus) + # The assessment in memory supplies what its saved file does, and every registered plot draws from both. + ef = res.error_fields + from_result, from_file = AEF._load(ef), AEF._load(h5path) + for k in keys(from_result) + k === :tolerances && continue + a, b = from_result[k], from_file[k] + @test (a === nothing && b === nothing) || (a isa AbstractArray{<:Number} ? isapprox(a, b; rtol=1e-10) : isequal(a, b)) + end + for e in GPEC.Analysis.PLOTS + GPEC.Analysis._accepts(e, ef) || continue + @test e.fn(ef) isa Plots.Plot + @test e.fn(h5path) isa Plots.Plot + end + # An in-memory table plots the same way a stored run does, which is the whole point of # the dual entry points: coil geometry that was never part of a run has no file to read. in_memory = AEF.plot_coil_sensitivities(table)