Analysis - BUGFIX - Read Jb_theta and Jb_zeta and smoke-test the gpec.h5 plots

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"]))

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"])

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")

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