From cd0adf64dba0bdec47251ae0376904aec25dea60 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 23 Oct 2026 20:34:26 +0200 Subject: [PATCH] Analysis - BUGFIX - Read Jb_theta, Jb_zeta and the kinetic rational_psi, and smoke-test the gpec.h5 plots --- src/Analysis/ForceFreeStates.jl | 2 +- src/Analysis/PerturbedEquilibrium.jl | 17 +++--- test/runtests.jl | 1 + test/runtests_analysis_plots.jl | 79 ++++++++++++++++++++++++++++ 4 files changed, 90 insertions(+), 9 deletions(-) create mode 100644 test/runtests_analysis_plots.jl diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index 304b9ba54..bff4d1dbd 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -465,7 +465,7 @@ function plot_cond_fbar(h5path; save_path=nothing, zoom=false) (read(kg["scan_psi"]), read(kg["scan_cond"]), read(kg["scan_threshold"]), - read(kg["psi"]), + read(kg["rational_psi"]), read(fid["SingularSurfaces/rational_psi"]), read(fid["SingularSurfaces/rational_q"]), read(kg["rational_count"])) diff --git a/src/Analysis/PerturbedEquilibrium.jl b/src/Analysis/PerturbedEquilibrium.jl index 33f8bdc0f..03f269660 100644 --- a/src/Analysis/PerturbedEquilibrium.jl +++ b/src/Analysis/PerturbedEquilibrium.jl @@ -342,8 +342,9 @@ Inspired by `plot_spectrograms.py` from OMFIT GPEC. ### Keyword arguments - `component`: Response field component to plot; one of `:xi_psi`, `:b_psi`, - `:b_theta`, `:b_zeta` (default: `:xi_psi`). `:b_psi` reads the area-normalized - `b_psi_area_weighted` dataset. + `: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 + contravariant field components. - `save_path`: If provided, save the figure to this path (default: `nothing`) ### Returns @@ -352,13 +353,13 @@ 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^ψ - :b_theta => "b_theta", - :b_zeta => "b_zeta", + :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, :b_theta, :b_zeta") + error("component must be one of :xi_psi, :b_psi, :Jb_theta, :Jb_zeta") base = "PerturbedEquilibrium/Response/" dataset_path = base * comp_map[component] @@ -368,7 +369,7 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) data, psi_response, mlow, mhigh, nhigh, q95, rational_psi = h5open(h5path, "r") do fid read(fid[dataset_path]), - read(fid["ForceFreeStates/Solutions/ForwardIntegration/psi"]), + 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"]) diff --git a/test/runtests.jl b/test/runtests.jl index b3891e63c..bc068435a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -50,6 +50,7 @@ else include("./runtests_kinetic_profiles.jl") include("./runtests_resist_eval.jl") include("./runtests_localstability.jl") + include("./runtests_analysis_plots.jl") include("./runtests_slayer_params.jl") include("./runtests_slayer_riccati.jl") include("./runtests_slayer_inputs.jl") diff --git a/test/runtests_analysis_plots.jl b/test/runtests_analysis_plots.jl new file mode 100644 index 000000000..63bd4b598 --- /dev/null +++ b/test/runtests_analysis_plots.jl @@ -0,0 +1,79 @@ +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 + 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}()) + template = joinpath(@__DIR__, "test_data", fixture) + for name in readdir(template) + cp(joinpath(template, name), joinpath(dir, name)) + end + toml_path = joinpath(dir, "gpec.toml") + inputs = TOML.parsefile(toml_path) + inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = true + merge!(inputs, sections) + open(io -> TOML.print(io, inputs), toml_path, "w") + GPEC.main([dir]) + return 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 + mktempdir() do dir + h5 = _run_fixture(dir, "regression_solovev_ideal_example"; sections=Dict{String,Any}( + "ForcingTerms" => Dict{String,Any}("forcing_data_format" => "coil", "mtheta_coil" => 240, "nzeta_coil" => 32, + "coil_set" => [Dict{String,Any}("name" => "hoop", "source" => "pf_hoop", "radius" => 1.5, "height" => 0.4, + "currents" => [2.0e3], "tiltx" => [3.0])]), + "PerturbedEquilibrium" => Dict{String,Any}("compute_response" => true, "compute_singular_coupling" => true, + "verbose" => false, "write_outputs_to_HDF5" => 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)) + end + for component in (:xi_psi, :b_psi, :Jb_theta, :Jb_zeta) + @test _drew(A.PerturbedEquilibrium.plot_mode_spectrogram(h5; component=component)) + end + @test_throws ErrorException 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 + + png = joinpath(dir, "qprofile.png") + A.Equilibrium.plot_qprofile(h5; save_path=png) + @test isfile(png) + + theta_data, theta_grid, n_vals = A.PerturbedEquilibriumModes.modes_to_theta(h5, "PerturbedEquilibrium/Response/xi_psi") + npsi = h5open(f -> length(read(f, "PerturbedEquilibrium/Response/psi")), h5) + @test size(theta_data) == (npsi, length(theta_grid), length(n_vals)) + end + end + + @testset "kinetic run draws the F̄ condition-number scan" begin + mktempdir() do dir + h5 = _run_fixture(dir, "regression_solovev_kinetic_example") + @test _drew(A.ForceFreeStates.plot_cond_fbar(h5)) + end + end +end