From b5f71a15aa48a6ec85334b281df4dedd6d38798f Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Sun, 22 Mar 2026 11:26:15 -0400 Subject: [PATCH 01/12] ANALYSIS - NEW FEATURE - Add equilibrium, stability, and PE visualization toolkit MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New component and summary plot functions across three Analysis submodules: ForceFreeStates (src/Analysis/ForceFreeStates.jl): - plot_energy_eigenvectors: heatmap of |Wt| eigenvectors vs (m, mode index) - plot_eigenvalue_spectrum: eigenvalue scatter colored by stability (red/green) - plot_delta_prime: |Δ'| bar chart from FFS asymptotic coefficients (ca_left/ca_right) - plot_dcon_summary: 4-panel composer (eigenvectors, criterion, spectrum, Δ') - plot_singular_surfaces: 2-panel composer (q-profile with surface markers, Δ') Equilibrium (src/Analysis/Equilibrium.jl): - plot_qprofile: q(ψ) with rational surface overlays - plot_pressure_profile: μ₀p(ψ) - plot_f_profile: 2πF(ψ) - plot_equilibrium_summary: 3-panel composer (q, pressure, F profiles) PerturbedEquilibrium (src/Analysis/PerturbedEquilibrium.jl, new file): - plot_resonant_flux, plot_island_widths, plot_chirikov_parameter, plot_pe_delta_prime - plot_resonant_field: 5-panel composer of resonant coupling quantities - plot_mode_spectrogram: 2-panel line+heatmap spectrogram (xi_psi/b_psi/b_theta/b_zeta) - plot_perturbed_equilibrium_summary: 3-panel composer (islands, energies, spectrogram) PE HDF5 output (src/PerturbedEquilibrium/Utils.jl): - Write resonant_flux, resonant_current, island_width_sq, penetrated_field, delta_prime (ComplexF64), island_half_width and chirikov_parameter (Float64) to singular_coupling group Updated all analyze_example.jl scripts to demonstrate the new functions. Co-Authored-By: Claude Sonnet 4.6 --- .../analyze_example.jl | 10 + .../Solovev_ideal_example/analyze_example.jl | 5 + .../analyze_example.jl | 5 + src/Analysis/Analysis.jl | 6 +- src/Analysis/Equilibrium.jl | 167 ++++++- src/Analysis/ForceFreeStates.jl | 266 +++++++++++ src/Analysis/PerturbedEquilibrium.jl | 425 ++++++++++++++++++ src/PerturbedEquilibrium/Utils.jl | 25 +- 8 files changed, 901 insertions(+), 8 deletions(-) create mode 100644 src/Analysis/PerturbedEquilibrium.jl diff --git a/examples/DIIID-like_ideal_example/analyze_example.jl b/examples/DIIID-like_ideal_example/analyze_example.jl index 192537588..7a562499b 100644 --- a/examples/DIIID-like_ideal_example/analyze_example.jl +++ b/examples/DIIID-like_ideal_example/analyze_example.jl @@ -8,3 +8,13 @@ h5path = joinpath(@__DIR__, "gpec.h5") p_modes = Analysis.ForceFreeStates.plot_mode_displacement(h5path; modes=1:5) p_eigen = Analysis.ForceFreeStates.plot_eigenmode_summary(h5path) p_stab = Analysis.ForceFreeStates.plot_stability_criterion(h5path) + +# Summary plots +p_dcon_summary = Analysis.ForceFreeStates.plot_dcon_summary(h5path) +p_sing = Analysis.ForceFreeStates.plot_singular_surfaces(h5path) +p_equil = Analysis.Equilibrium.plot_equilibrium_summary(h5path) + +# Perturbed equilibrium plots (requires PE module run with forcing.dat) +p_resfield = Analysis.PerturbedEquilibrium.plot_resonant_field(h5path) +p_pe_summary = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path) +p_spectrogram = Analysis.PerturbedEquilibrium.plot_mode_spectrogram(h5path; component=:xi_psi) diff --git a/examples/Solovev_ideal_example/analyze_example.jl b/examples/Solovev_ideal_example/analyze_example.jl index 192537588..45a2754e4 100644 --- a/examples/Solovev_ideal_example/analyze_example.jl +++ b/examples/Solovev_ideal_example/analyze_example.jl @@ -8,3 +8,8 @@ h5path = joinpath(@__DIR__, "gpec.h5") p_modes = Analysis.ForceFreeStates.plot_mode_displacement(h5path; modes=1:5) p_eigen = Analysis.ForceFreeStates.plot_eigenmode_summary(h5path) p_stab = Analysis.ForceFreeStates.plot_stability_criterion(h5path) + +# Summary plots +p_dcon_summary = Analysis.ForceFreeStates.plot_dcon_summary(h5path) +p_sing = Analysis.ForceFreeStates.plot_singular_surfaces(h5path) +p_equil = Analysis.Equilibrium.plot_equilibrium_summary(h5path) diff --git a/examples/Solovev_ideal_example_multi_n/analyze_example.jl b/examples/Solovev_ideal_example_multi_n/analyze_example.jl index 724ab4d35..b4e4e2539 100644 --- a/examples/Solovev_ideal_example_multi_n/analyze_example.jl +++ b/examples/Solovev_ideal_example_multi_n/analyze_example.jl @@ -10,6 +10,11 @@ p_modes = Analysis.ForceFreeStates.plot_mode_displacement(h5path; modes=1:5) p_eigen = Analysis.ForceFreeStates.plot_eigenmode_summary(h5path) p_stab = Analysis.ForceFreeStates.plot_stability_criterion(h5path) +# Summary plots +p_dcon_summary = Analysis.ForceFreeStates.plot_dcon_summary(h5path) +p_sing = Analysis.ForceFreeStates.plot_singular_surfaces(h5path) +p_equil = Analysis.Equilibrium.plot_equilibrium_summary(h5path) + # Analyze the single-n runs h5path_n1 = joinpath(@__DIR__, "single_n_1", "euler_n1.h5") h5path_n2 = joinpath(@__DIR__, "single_n_2", "euler_n2.h5") diff --git a/src/Analysis/Analysis.jl b/src/Analysis/Analysis.jl index 1a770894e..0a2a5765f 100644 --- a/src/Analysis/Analysis.jl +++ b/src/Analysis/Analysis.jl @@ -5,12 +5,14 @@ Post-processing and visualization utilities for GPEC simulation outputs. ## Submodules - - `ForceFreeStates`: Plotting functions for DCON-style ideal MHD stability results - - `Equilibrium`: Plotting functions for equilibrium objects + - `ForceFreeStates`: Plotting functions for ForceFreeStates (DCON-style ideal MHD stability) results + - `Equilibrium`: Plotting functions for equilibrium profiles and flux surfaces + - `PerturbedEquilibrium`: Plotting functions for perturbed equilibrium and singular coupling results """ module Analysis include("ForceFreeStates.jl") include("Equilibrium.jl") +include("PerturbedEquilibrium.jl") end # module Analysis diff --git a/src/Analysis/Equilibrium.jl b/src/Analysis/Equilibrium.jl index 19e745c8c..12380dae8 100644 --- a/src/Analysis/Equilibrium.jl +++ b/src/Analysis/Equilibrium.jl @@ -1,10 +1,11 @@ """ Equilibrium -Post-processing and visualization functions for GPEC equilibrium objects. +Post-processing and visualization functions for GPEC equilibrium objects and HDF5 outputs. """ module Equilibrium +using HDF5 using Plots """ @@ -66,4 +67,168 @@ function plot_flux_surfaces(plasma_eq; n_psi=11, n_theta=13) return p end +""" + plot_qprofile(h5path; 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 + +### Keyword arguments + + - `show_singular`: If `true`, overlay rational surface locations (default: `true`) + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +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["splines/profiles/xs"]), read(fid["splines/profiles/q"]), + read(fid["equil/q0"]), read(fid["equil/q95"]) + end + + p = plot( + xs, q; + xlabel="ψ_N", + ylabel="q", + title="Safety factor q(ψ)", + legend=false + ) + hline!(p, [q0, q95]; linestyle=:dot, color=:gray, label=nothing) + annotate!(p, 1.0, q0, text(" q0=$(round(q0, digits=2))", 7, :left, :gray)) + 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["singular/msing"]), read(fid["singular/psi"]), read(fid["singular/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)) + end + end + + isnothing(save_path) || savefig(p, save_path) + return p +end + +""" + plot_pressure_profile(h5path; 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 + +### Keyword arguments + + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_pressure_profile(h5path; save_path=nothing) + xs, mu0p, betat, betap1, msing, psi_sing = h5open(h5path, "r") do fid + read(fid["splines/profiles/xs"]), read(fid["splines/profiles/mu0p"]), + read(fid["equil/betat"]), read(fid["equil/betap1"]), + read(fid["singular/msing"]), read(fid["singular/psi"]) + end + + p = plot( + xs, mu0p; + xlabel="ψ_N", + ylabel="μ₀p", + title="Pressure profile (βₜ=$(round(betat, digits=3)), βₚ=$(round(betap1, digits=3)))", + legend=false + ) + for s in 1:msing + vline!(p, [psi_sing[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 the toroidal field function 2πF(ψ) profile (F = RBφ/(2π)). + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file + +### Keyword arguments + + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_f_profile(h5path; save_path=nothing) + xs, twopif = h5open(h5path, "r") do fid + read(fid["splines/profiles/xs"]), read(fid["splines/profiles/2piF"]) + end + + p = plot( + xs, twopif; + xlabel="ψ_N", + ylabel="2πF", + title="Toroidal field function 2πF(ψ)", + legend=false + ) + + isnothing(save_path) || savefig(p, save_path) + return p +end + +""" + plot_equilibrium_summary(h5path; save_path=nothing) + +Three-panel summary of equilibrium profiles: + + - q(ψ) safety factor with rational surface markers (`plot_qprofile`) + - μ₀p(ψ) pressure profile (`plot_pressure_profile`) + - 2πF(ψ) toroidal field function (`plot_f_profile`) + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file + +### Keyword arguments + + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +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["equil/q0"]), read(fid["equil/q95"]), + read(fid["equil/betat"]), read(fid["equil/betan"]), + read(fid["equil/kappa"]), read(fid["equil/li1"]) + end + + title_str = "q0=$(round(q0,digits=2)) q95=$(round(q95,digits=2)) βₜ=$(round(betat,digits=3)) βₙ=$(round(betan,digits=3)) κ=$(round(kappa,digits=2)) li1=$(round(li1,digits=3))" + + p_q = plot_qprofile(h5path; show_singular=true) + p_pres = plot_pressure_profile(h5path) + p_f = plot_f_profile(h5path) + + p = plot(p_q, p_pres, p_f; layout=(1, 3), size=(1200, 400), plot_title=title_str) + + isnothing(save_path) || savefig(p, save_path) + return p +end + end # module Equilibrium diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index 1fde16391..749d7b2bb 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -156,4 +156,270 @@ function plot_stability_criterion(h5path; save_path=nothing) return p end +""" + plot_energy_eigenvectors(h5path; 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 +`vacuum/wt`). Plasma and vacuum eigenvectors are not stored separately in the HDF5 output. + +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`) + +### Keyword arguments + + - `matrix_type`: Energy matrix to plot; only `:total` is currently supported + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_energy_eigenvectors(h5path; 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["vacuum/wt"]), read(fid["equil/psio"]), read(fid["info/mlow"]) + end + + isempty(wt) && error("No vacuum data in $h5path; rerun with vac_flag = true") + + chi1 = 2π * psio + wt_scaled = wt * (chi1 * 1e-3) + + nmn = size(wt_scaled, 1) + nmodes = size(wt_scaled, 2) + m_vals = (0:(nmn-1)) .+ mlow + + p = heatmap( + m_vals, 1:nmodes, abs.(wt_scaled'); + xlabel="m", + ylabel="mode index", + title="Energy eigenvectors |Wₜ| (total)", + colorbar_title="|Wₜ|" + ) + + isnothing(save_path) || savefig(p, save_path) + return p +end + +""" + plot_eigenvalue_spectrum(h5path; 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. + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file with vacuum data (`vac_flag = true`) + +### Keyword arguments + + - `matrix_type`: Which eigenvalues to plot: `:total` (`et`), `:plasma` (`ep`), or `:vacuum` (`ev`) + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_eigenvalue_spectrum(h5path; matrix_type=:total, save_path=nothing) + dataset = Dict(:total => "vacuum/et", :plasma => "vacuum/ep", :vacuum => "vacuum/ev") + 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") + + nmodes = length(et) + ev_real = real.(et) + colors = [v > 0 ? :red : :green for v in ev_real] + + p = scatter( + ev_real, 1:nmodes; + xlabel="Re(eigenvalue)", + ylabel="mode index", + title="Eigenvalue spectrum ($matrix_type)", + legend=false, + color=colors, + markerstrokewidth=0 + ) + hline!(p, [0]; linestyle=:dash, color=:black, label=nothing) + + isnothing(save_path) || savefig(p, save_path) + return p +end + +""" + plot_delta_prime(h5path; save_path=nothing) + +Bar chart of `|Δ'|` per singular surface, computed from the stored asymptotic coefficients +`ca_left` and `ca_right`. Bars are colored red (tearing unstable, Re(Δ') > 0) or blue +(tearing stable). + +Δ' 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`. + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file + +### Keyword arguments + + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_delta_prime(h5path; save_path=nothing) + msing, q_sing, ca_l, ca_r, psio, mn_index = h5open(h5path, "r") do fid + read(fid["singular/msing"]), read(fid["singular/q"]), + read(fid["singular/ca_left"]), read(fid["singular/ca_right"]), + read(fid["equil/psio"]), read(fid["info/mn_index"]) + end + + msing == 0 && return plot(; title="No singular surfaces found", legend=false) + + numpert_total = size(ca_l, 1) + chi1 = 2π * psio + + dp_vals = ComplexF64[] + labels = String[] + + for s in 1:msing + q_s = q_sing[s] + # Find the resonant mode index: (m, n) in mn_index with m/n ≈ q_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 + push!(labels, "q=$(round(q_s, digits=3))") + end + + colors = [real(v) > 0 ? :red : :steelblue for v in dp_vals] + + p = bar( + 1:msing, abs.(dp_vals); + xticks=(1:msing, labels), + xlabel="rational surface", + ylabel="|Δ'|", + title="Tearing stability Δ' (FFS asymptotic coefficients)", + legend=false, + color=colors, + xrotation=30 + ) + + isnothing(save_path) || savefig(p, save_path) + return p +end + +""" + plot_dcon_summary(h5path; save_path=nothing) + +Four-panel summary of ForceFreeStates (DCON-style) stability results, combining: + + - Energy eigenvector heatmap (`plot_energy_eigenvectors`) + - Stability criterion vs ψ_N (`plot_stability_criterion`) + - Eigenvalue spectrum (`plot_eigenvalue_spectrum`) + - 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. + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file + +### Keyword arguments + + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_dcon_summary(h5path; save_path=nothing) + has_vac = h5open(h5path, "r") do fid + haskey(fid, "vacuum/wt") && !isempty(read(fid["vacuum/wt"])) + end + + p_crit = plot_stability_criterion(h5path) + p_dp = plot_delta_prime(h5path) + + if has_vac + p_evec = plot_energy_eigenvectors(h5path; matrix_type=:total) + p_eval = plot_eigenvalue_spectrum(h5path; matrix_type=:total) + p = plot(p_evec, p_crit, p_eval, p_dp; layout=(2, 2), size=(1000, 800)) + else + annotate!(p_crit, :center, :top, + text("No vacuum data — rerun with vac_flag = true for full summary", 8, :gray)) + p = plot(p_crit, p_dp; layout=(1, 2), size=(1000, 400)) + end + + isnothing(save_path) || savefig(p, save_path) + return p +end + +""" + plot_singular_surfaces(h5path; save_path=nothing) + +Two-panel summary of singular (rational) surface locations and tearing stability: + + - Left: q(ψ) profile with vertical markers at each rational surface + - Right: `|Δ'|` bar chart at each rational surface (`plot_delta_prime`) + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file + +### Keyword arguments + + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_singular_surfaces(h5path; save_path=nothing) + xs, q_profile, q0, q95, msing, psi_sing, q_sing = h5open(h5path, "r") do fid + read(fid["splines/profiles/xs"]), read(fid["splines/profiles/q"]), + read(fid["equil/q0"]), read(fid["equil/q95"]), + read(fid["singular/msing"]), read(fid["singular/psi"]), read(fid["singular/q"]) + end + + p_q = plot( + xs, q_profile; + xlabel="ψ_N", + ylabel="q", + title="Safety factor and rational surfaces", + legend=false + ) + hline!(p_q, [q0, q95]; linestyle=:dot, color=:gray, label=nothing) + for s in 1:msing + vline!(p_q, [psi_sing[s]]; linestyle=:dash, color=:red, label=nothing) + annotate!(p_q, psi_sing[s], q_sing[s], + text(" q=$(round(q_sing[s], digits=2))", 7, :left, :red)) + end + + p_dp = plot_delta_prime(h5path) + + p = plot(p_q, p_dp; layout=(1, 2), size=(1000, 400)) + + isnothing(save_path) || savefig(p, save_path) + return p +end + end # module ForceFreeStates diff --git a/src/Analysis/PerturbedEquilibrium.jl b/src/Analysis/PerturbedEquilibrium.jl new file mode 100644 index 000000000..f261acfbd --- /dev/null +++ b/src/Analysis/PerturbedEquilibrium.jl @@ -0,0 +1,425 @@ +""" + PerturbedEquilibrium + +Post-processing and visualization functions for GPEC perturbed equilibrium results stored +in the `perturbed_equilibrium/` group of a GPEC HDF5 output file. +""" +module PerturbedEquilibrium + +using HDF5 +using Plots + +# Check that a PE dataset exists and is non-empty. +function _has_pe_data(h5path, key) + h5open(h5path, "r") do fid + haskey(fid, key) && !isempty(read(fid[key])) + end +end + +""" + plot_resonant_flux(h5path; save_path=nothing) + +Bar chart of `|Φ_res|` (normalized resonant flux) per singular surface. One bar series per +toroidal mode n, labeled by n value. + +Requires the perturbed equilibrium module to have been run and `singular_coupling/resonant_flux` +to be present in the HDF5 file. + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file with perturbed equilibrium output + +### Keyword arguments + + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_resonant_flux(h5path; save_path=nothing) + key = "perturbed_equilibrium/singular_coupling/resonant_flux" + _has_pe_data(h5path, key) || + return plot(; title="No resonant flux data — run with perturbed equilibrium enabled", legend=false) + + resonant_flux, q_sing, msing, mn_index = h5open(h5path, "r") do fid + read(fid[key]), read(fid["singular/q"]), + read(fid["singular/msing"]), read(fid["info/mn_index"]) + end + + labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] + xticks_arg = (1:msing, labels) + + p = plot(; xlabel="rational surface", ylabel="|Φ_res|", + title="Resonant flux |Φ_res| per surface", legend=:outertopright) + + n_vals = unique(mn_index[:, 2]) + for nn in n_vals + n_rows = findall(j -> mn_index[j, 2] == nn, 1:size(mn_index, 1)) + # Sum over poloidal modes for this n + rf_n = [sum(abs.(resonant_flux[n_rows, s])) for s in 1:msing] + bar!(p, 1:msing, rf_n; label="n=$nn", alpha=0.7) + end + plot!(p; xticks=xticks_arg, xrotation=30) + + isnothing(save_path) || savefig(p, save_path) + return p +end + +""" + plot_island_widths(h5path; save_path=nothing) + +Bar chart of island half-width `w/2` per singular surface. + +Requires `singular_coupling/island_half_width` in the HDF5 file. + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file with perturbed equilibrium output + +### Keyword arguments + + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_island_widths(h5path; save_path=nothing) + key = "perturbed_equilibrium/singular_coupling/island_half_width" + _has_pe_data(h5path, key) || + return plot(; title="No island width data — run with perturbed equilibrium enabled", legend=false) + + island_hw, q_sing, msing = h5open(h5path, "r") do fid + read(fid[key]), read(fid["singular/q"]), read(fid["singular/msing"]) + end + + labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] + + p = bar( + 1:msing, island_hw; + xticks=(1:msing, labels), + xlabel="rational surface", + ylabel="w/2", + title="Island half-widths", + legend=false, + color=:steelblue, + xrotation=30 + ) + + isnothing(save_path) || savefig(p, save_path) + return p +end + +""" + plot_chirikov_parameter(h5path; save_path=nothing) + +Bar chart of the Chirikov overlap parameter per singular surface, with a horizontal reference +line at K = 1 (island overlap threshold). Bars are colored red when K > 1. + +Requires `singular_coupling/chirikov_parameter` in the HDF5 file. + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file with perturbed equilibrium output + +### Keyword arguments + + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_chirikov_parameter(h5path; save_path=nothing) + key = "perturbed_equilibrium/singular_coupling/chirikov_parameter" + _has_pe_data(h5path, key) || + return plot(; title="No Chirikov data — run with perturbed equilibrium enabled", legend=false) + + chirikov, q_sing, msing = h5open(h5path, "r") do fid + read(fid[key]), read(fid["singular/q"]), read(fid["singular/msing"]) + end + + labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] + colors = [k > 1.0 ? :red : :steelblue for k in chirikov] + + p = bar( + 1:msing, chirikov; + xticks=(1:msing, labels), + xlabel="rational surface", + ylabel="K (Chirikov)", + title="Chirikov overlap parameter", + legend=false, + color=colors, + xrotation=30 + ) + hline!(p, [1.0]; linestyle=:dash, color=:black, label=nothing) + + isnothing(save_path) || savefig(p, save_path) + return p +end + +""" + plot_pe_delta_prime(h5path; save_path=nothing) + +Bar chart of `|Δ'|` per singular surface computed by the perturbed equilibrium module +(from `singular_coupling/delta_prime`). One bar series per toroidal mode n. + +This is complementary to `Analysis.ForceFreeStates.plot_delta_prime`, which uses the FFS +asymptotic coefficients. The PE result includes the vacuum Green's function contribution. + +Requires `singular_coupling/delta_prime` in the HDF5 file. + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file with perturbed equilibrium output + +### Keyword arguments + + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_pe_delta_prime(h5path; save_path=nothing) + key = "perturbed_equilibrium/singular_coupling/delta_prime" + _has_pe_data(h5path, key) || + return plot(; title="No PE Δ' data — run with perturbed equilibrium enabled", legend=false) + + delta_prime, q_sing, msing, mn_index = h5open(h5path, "r") do fid + read(fid[key]), read(fid["singular/q"]), read(fid["singular/msing"]), + read(fid["info/mn_index"]) + end + + labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] + + p = plot(; xlabel="rational surface", ylabel="|Δ'|", + title="Tearing stability Δ' (PE)", legend=:outertopright) + + n_vals = unique(mn_index[:, 2]) + for nn in n_vals + n_rows = findall(j -> mn_index[j, 2] == nn, 1:size(mn_index, 1)) + dp_n = [maximum(abs.(delta_prime[n_rows, s])) for s in 1:msing] + bar!(p, 1:msing, dp_n; label="n=$nn", alpha=0.7) + end + plot!(p; xticks=(1:msing, labels), xrotation=30) + + isnothing(save_path) || savefig(p, save_path) + return p +end + +""" + plot_resonant_field(h5path; save_path=nothing) + +Five-panel summary of resonant coupling quantities at each singular surface: + + - `|Φ_res|`: resonant flux (`plot_resonant_flux`) + - `|Δ'|`: tearing stability parameter (`plot_pe_delta_prime`) + - `|I_res|`: resonant current + - `w/2`: island half-width (`plot_island_widths`) + - `K`: Chirikov overlap parameter (`plot_chirikov_parameter`) + +Inspired by `plot_resonant_field.py` from OMFIT GPEC. + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file with perturbed equilibrium output + +### Keyword arguments + + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_resonant_field(h5path; save_path=nothing) + p1 = plot_resonant_flux(h5path) + p2 = plot_pe_delta_prime(h5path) + p3 = _plot_resonant_current(h5path) + p4 = plot_island_widths(h5path) + p5 = plot_chirikov_parameter(h5path) + + p = plot(p1, p2, p3, p4, p5; layout=(5, 1), size=(800, 1400)) + + isnothing(save_path) || savefig(p, save_path) + return p +end + +# Internal helper — resonant current bar chart +function _plot_resonant_current(h5path) + key = "perturbed_equilibrium/singular_coupling/resonant_current" + _has_pe_data(h5path, key) || + return plot(; title="No resonant current data", legend=false) + + resonant_current, q_sing, msing, mn_index = h5open(h5path, "r") do fid + read(fid[key]), read(fid["singular/q"]), read(fid["singular/msing"]), + read(fid["info/mn_index"]) + end + + labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] + + p = plot(; xlabel="rational surface", ylabel="|I_res|", + title="Resonant current |I_res| per surface", legend=:outertopright) + + n_vals = unique(mn_index[:, 2]) + for nn in n_vals + n_rows = findall(j -> mn_index[j, 2] == nn, 1:size(mn_index, 1)) + rc_n = [sum(abs.(resonant_current[n_rows, s])) for s in 1:msing] + bar!(p, 1:msing, rc_n; label="n=$nn", alpha=0.7) + end + plot!(p; xticks=(1:msing, labels), xrotation=30) + + return p +end + +""" + plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) + +Two-panel spectrogram of a perturbed equilibrium response field component: + + - Top: `|component|` vs ψ_N, one curve per poloidal mode m (colored by m) + - Bottom: Heatmap of `|component|` in (ψ_N × m) space, with white dashed lines at + rational surface locations + +Inspired by `plot_spectrograms.py` from OMFIT GPEC. + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file with perturbed equilibrium response output + +### Keyword arguments + + - `component`: Response field component to plot; one of `:xi_psi`, `:b_psi`, + `:b_theta`, `:b_zeta` (default: `:xi_psi`) + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) + comp_map = Dict( + :xi_psi => ("xi_psi_real", "xi_psi_imag"), + :b_psi => ("b_psi_real", "b_psi_imag"), + :b_theta => ("b_theta_real", "b_theta_imag"), + :b_zeta => ("b_zeta_real", "b_zeta_imag"), + ) + haskey(comp_map, component) || + error("component must be one of :xi_psi, :b_psi, :b_theta, :b_zeta") + + rkey, ikey = comp_map[component] + base = "perturbed_equilibrium/response/" + + _has_pe_data(h5path, base * rkey) || + return plot(; title="No response data — run with perturbed equilibrium enabled", legend=false) + + data_r, data_i, xs, mlow, mhigh, msing, psi_sing = h5open(h5path, "r") do fid + read(fid[base * rkey]), read(fid[base * ikey]), + read(fid["splines/profiles/xs"]), + read(fid["info/mlow"]), read(fid["info/mhigh"]), + read(fid["singular/msing"]), read(fid["singular/psi"]) + end + + data = complex.(data_r, data_i) # shape: (numpert_total, npsi) + mpert = mhigh - mlow + 1 + + # Use only the first n's modes for a clean spectrogram (single-n assumption for display) + m_vals = mlow:mhigh + data_mn = data[1:mpert, :] # first mpert rows correspond to first n + + # Top panel: line plot per mode + p1 = plot(; + xlabel="ψ_N", + ylabel="|$(component)|", + title="Mode spectrogram: $(component)", + legend=:outertopright + ) + cmap = cgrad(:roma, mpert; categorical=true) + for (i, m) in enumerate(m_vals) + plot!(p1, xs, abs.(data_mn[i, :]); label="m=$m", color=cmap[i], linewidth=1.5) + end + + # Bottom panel: heatmap + p2 = heatmap( + xs, collect(m_vals), abs.(data_mn); + xlabel="ψ_N", + ylabel="m", + title="", + colorbar_title="|$(component)|" + ) + # Overlay rational surface locations as white dashed lines + for s in 1:msing + vline!(p2, [psi_sing[s]]; linestyle=:dash, color=:white, linewidth=1.5, label=nothing) + end + + p = plot(p1, p2; layout=(2, 1), size=(900, 700)) + + isnothing(save_path) || savefig(p, save_path) + return p +end + +""" + plot_perturbed_equilibrium_summary(h5path; save_path=nothing) + +Three-panel composite summary of perturbed equilibrium results: + + - Top-left: Island half-widths and Chirikov parameter overlay (`plot_island_widths` + + `plot_chirikov_parameter` on shared axes) + - Top-right: Energy breakdown — plasma, vacuum, and total energies + - Bottom: ξ_ψ mode spectrogram (`plot_mode_spectrogram`) + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file with perturbed equilibrium output + +### Keyword arguments + + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_perturbed_equilibrium_summary(h5path; save_path=nothing) + p_islands = plot_island_widths(h5path) + p_energies = _plot_energies(h5path) + p_spectro = plot_mode_spectrogram(h5path; component=:xi_psi) + + l = @layout [grid(1, 2){0.35h}; b] + p = plot(p_islands, p_energies, p_spectro; layout=l, size=(1000, 1000)) + + isnothing(save_path) || savefig(p, save_path) + return p +end + +# Internal helper — energy breakdown bar chart +function _plot_energies(h5path) + key = "perturbed_equilibrium/energies/plasma_energy" + _has_pe_data(h5path, key) || + return plot(; title="No energy data", legend=false) + + ep, ev, et = h5open(h5path, "r") do fid + read(fid["perturbed_equilibrium/energies/plasma_energy"]), + read(fid["perturbed_equilibrium/energies/vacuum_energy"]), + read(fid["perturbed_equilibrium/energies/total_energy"]) + end + + vals = [real(ep), real(ev), real(et)] + names = ["Plasma", "Vacuum", "Total"] + colors = [:steelblue, :darkorange, :green] + + p = bar( + names, vals; + ylabel="Energy", + title="Energy breakdown", + legend=false, + color=colors + ) + hline!(p, [0]; linestyle=:dash, color=:black, label=nothing) + + return p +end + +end # module PerturbedEquilibrium diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl index 03aea908d..652f1023b 100644 --- a/src/PerturbedEquilibrium/Utils.jl +++ b/src/PerturbedEquilibrium/Utils.jl @@ -61,12 +61,20 @@ perturbed_equilibrium/ │ ├── amplitude_real # Real parts of forcing amplitudes │ └── amplitude_imag # Imaginary parts of forcing amplitudes ├── response/ -│ ├── xi_perturbed # Displacement field -│ └── b_perturbed # Magnetic field perturbation +│ ├── xi_psi_real/imag # Radial displacement (real/imag parts) +│ ├── b_psi_real/imag # Normal field component +│ ├── b_theta_real/imag # Poloidal field component +│ └── b_zeta_real/imag # Toroidal field component ├── singular_coupling/ -│ ├── coupling_coefficient_real -│ ├── coupling_coefficient_imag -│ └── resonant_amplitude +│ ├── coupling_coefficient_real/imag +│ ├── resonant_amplitude +│ ├── resonant_flux # ComplexF64 [numpert_total × msing] +│ ├── resonant_current # ComplexF64 [numpert_total × msing] +│ ├── island_width_sq # ComplexF64 [numpert_total × msing] +│ ├── penetrated_field # ComplexF64 [numpert_total × msing] +│ ├── delta_prime # ComplexF64 [numpert_total × msing]; tearing stability Δ' +│ ├── island_half_width # Float64 [msing]; actual w/2 +│ └── chirikov_parameter # Float64 [msing]; island overlap metric └── energies/ ├── plasma_energy ├── vacuum_energy @@ -113,6 +121,13 @@ function write_outputs_to_HDF5( coupling_group["coupling_coefficient_real"] = real(state.coupling_coefficient) coupling_group["coupling_coefficient_imag"] = imag(state.coupling_coefficient) coupling_group["resonant_amplitude"] = state.resonant_amplitude + coupling_group["resonant_flux"] = state.resonant_flux + coupling_group["resonant_current"] = state.resonant_current + coupling_group["island_width_sq"] = state.island_width_sq + coupling_group["penetrated_field"] = state.penetrated_field + coupling_group["delta_prime"] = state.delta_prime + coupling_group["island_half_width"] = state.island_half_width + coupling_group["chirikov_parameter"] = state.chirikov_parameter # Write additional metrics for (key, val) in intr.singular_coupling_metrics From bce236833a5ae4af1df4f12181559fbf0ad34d7e Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Mon, 23 Mar 2026 10:24:34 -0400 Subject: [PATCH 02/12] ANALYSIS - BUG FIX - Fix plot bugs: vline/hline, annotate, PE mn_index, response shape - ForceFreeStates: Replace hline! with vline! in plot_eigenvalue_spectrum (stability boundary is Re(eigenvalue)=0, a vertical line, not horizontal) - ForceFreeStates: Replace invalid annotate!(p, :center, :top, ...) with title! in plot_dcon_summary (no-vacuum-data branch); :center/:top are not data coords - ForceFreeStates: Add inline note explaining imag(xi_psi) phase convention - PerturbedEquilibrium: Fix plot_resonant_flux, plot_pe_delta_prime, _plot_resonant_current to read PE forcing mode n from perturbed_equilibrium/forcing_modes/n instead of info/mn_index (FFS). PE coupling matrices are (num_forcing_modes, msing), not (numpert_total_FFS, msing) - PerturbedEquilibrium: Fix plot_mode_spectrogram data layout and psi axis. Response data is (npsi, numpert_total) not (numpert_total, npsi); psi axis is integration/psi (737 pts), not splines/profiles/xs (183 pts) - PerturbedEquilibrium: Fix plot_perturbed_equilibrium_summary docstring to match implementation (no Chirikov overlay on island-width panel) - Add examples/DIIID-like_ideal_example/make_plots.jl: saves all 20 plots to examples/DIIID-like_ideal_example/plots/ using GR backend Co-Authored-By: Claude Sonnet 4.6 --- .../DIIID-like_ideal_example/make_plots.jl | 73 +++++++++++++++++++ src/Analysis/ForceFreeStates.jl | 7 +- src/Analysis/PerturbedEquilibrium.jl | 42 +++++------ 3 files changed, 97 insertions(+), 25 deletions(-) create mode 100644 examples/DIIID-like_ideal_example/make_plots.jl diff --git a/examples/DIIID-like_ideal_example/make_plots.jl b/examples/DIIID-like_ideal_example/make_plots.jl new file mode 100644 index 000000000..63c6f9f66 --- /dev/null +++ b/examples/DIIID-like_ideal_example/make_plots.jl @@ -0,0 +1,73 @@ +using Pkg +Pkg.activate(joinpath(@__DIR__, "../..")) +using GeneralizedPerturbedEquilibrium, Plots, GeneralizedPerturbedEquilibrium.Analysis +gr() + +h5path = joinpath(@__DIR__, "gpec.h5") +outdir = joinpath(@__DIR__, "plots") +mkpath(outdir) + +# --- Equilibrium --- +p = Analysis.Equilibrium.plot_equilibrium_summary(h5path) +savefig(p, joinpath(outdir, "equil_summary.png")) + +p = Analysis.Equilibrium.plot_qprofile(h5path) +savefig(p, joinpath(outdir, "equil_qprofile.png")) + +p = Analysis.Equilibrium.plot_pressure_profile(h5path) +savefig(p, joinpath(outdir, "equil_pressure.png")) + +p = Analysis.Equilibrium.plot_f_profile(h5path) +savefig(p, joinpath(outdir, "equil_F.png")) + +# --- ForceFreeStates --- +p = Analysis.ForceFreeStates.plot_mode_displacement(h5path; modes=1:5) +savefig(p, joinpath(outdir, "ffs_mode_displacement.png")) + +p = Analysis.ForceFreeStates.plot_eigenmode_summary(h5path) +savefig(p, joinpath(outdir, "ffs_eigenmode_summary.png")) + +p = Analysis.ForceFreeStates.plot_stability_criterion(h5path) +savefig(p, joinpath(outdir, "ffs_stability_criterion.png")) + +p = Analysis.ForceFreeStates.plot_energy_eigenvectors(h5path) +savefig(p, joinpath(outdir, "ffs_energy_eigenvectors.png")) + +p = Analysis.ForceFreeStates.plot_eigenvalue_spectrum(h5path) +savefig(p, joinpath(outdir, "ffs_eigenvalue_spectrum.png")) + +p = Analysis.ForceFreeStates.plot_delta_prime(h5path) +savefig(p, joinpath(outdir, "ffs_delta_prime.png")) + +p = Analysis.ForceFreeStates.plot_dcon_summary(h5path) +savefig(p, joinpath(outdir, "ffs_dcon_summary.png")) + +p = Analysis.ForceFreeStates.plot_singular_surfaces(h5path) +savefig(p, joinpath(outdir, "ffs_singular_surfaces.png")) + +# --- PerturbedEquilibrium --- +p = Analysis.PerturbedEquilibrium.plot_resonant_flux(h5path) +savefig(p, joinpath(outdir, "pe_resonant_flux.png")) + +p = Analysis.PerturbedEquilibrium.plot_island_widths(h5path) +savefig(p, joinpath(outdir, "pe_island_widths.png")) + +p = Analysis.PerturbedEquilibrium.plot_chirikov_parameter(h5path) +savefig(p, joinpath(outdir, "pe_chirikov.png")) + +p = Analysis.PerturbedEquilibrium.plot_pe_delta_prime(h5path) +savefig(p, joinpath(outdir, "pe_delta_prime.png")) + +p = Analysis.PerturbedEquilibrium.plot_resonant_field(h5path) +savefig(p, joinpath(outdir, "pe_resonant_field.png")) + +p = Analysis.PerturbedEquilibrium.plot_mode_spectrogram(h5path; component=:xi_psi) +savefig(p, joinpath(outdir, "pe_spectrogram_xi_psi.png")) + +p = Analysis.PerturbedEquilibrium.plot_mode_spectrogram(h5path; component=:b_psi) +savefig(p, joinpath(outdir, "pe_spectrogram_b_psi.png")) + +p = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path) +savefig(p, joinpath(outdir, "pe_summary.png")) + +println("All plots saved to: $outdir") diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index 749d7b2bb..89e53758b 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -43,7 +43,7 @@ function plot_mode_displacement(h5path; modes=1:5, save_path=nothing) ) for m in modes mlow <= m <= mhigh || continue - plot!(p, psi, imag.(xi_psi[m-mlow+1, 1, :]); label="m=$m") + plot!(p, psi, imag.(xi_psi[m-mlow+1, 1, :]); label="m=$m") # DCON phase convention: ξ_ψ is purely imaginary on the real axis end isnothing(save_path) || savefig(p, save_path) @@ -250,7 +250,7 @@ function plot_eigenvalue_spectrum(h5path; matrix_type=:total, save_path=nothing) color=colors, markerstrokewidth=0 ) - hline!(p, [0]; linestyle=:dash, color=:black, label=nothing) + vline!(p, [0]; linestyle=:dash, color=:black, label=nothing) isnothing(save_path) || savefig(p, save_path) return p @@ -364,8 +364,7 @@ function plot_dcon_summary(h5path; save_path=nothing) p_eval = plot_eigenvalue_spectrum(h5path; matrix_type=:total) p = plot(p_evec, p_crit, p_eval, p_dp; layout=(2, 2), size=(1000, 800)) else - annotate!(p_crit, :center, :top, - text("No vacuum data — rerun with vac_flag = true for full summary", 8, :gray)) + title!(p_crit, "Stability criterion (no vacuum data — rerun with vac_flag = true)") p = plot(p_crit, p_dp; layout=(1, 2), size=(1000, 400)) end diff --git a/src/Analysis/PerturbedEquilibrium.jl b/src/Analysis/PerturbedEquilibrium.jl index f261acfbd..f77a06d5e 100644 --- a/src/Analysis/PerturbedEquilibrium.jl +++ b/src/Analysis/PerturbedEquilibrium.jl @@ -42,9 +42,10 @@ function plot_resonant_flux(h5path; save_path=nothing) _has_pe_data(h5path, key) || return plot(; title="No resonant flux data — run with perturbed equilibrium enabled", legend=false) - resonant_flux, q_sing, msing, mn_index = h5open(h5path, "r") do fid + resonant_flux, q_sing, msing, pe_n = h5open(h5path, "r") do fid read(fid[key]), read(fid["singular/q"]), - read(fid["singular/msing"]), read(fid["info/mn_index"]) + read(fid["singular/msing"]), + read(fid["perturbed_equilibrium/forcing_modes/n"]) end labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] @@ -53,9 +54,9 @@ function plot_resonant_flux(h5path; save_path=nothing) p = plot(; xlabel="rational surface", ylabel="|Φ_res|", title="Resonant flux |Φ_res| per surface", legend=:outertopright) - n_vals = unique(mn_index[:, 2]) + n_vals = unique(pe_n) for nn in n_vals - n_rows = findall(j -> mn_index[j, 2] == nn, 1:size(mn_index, 1)) + n_rows = findall(==(nn), pe_n) # Sum over poloidal modes for this n rf_n = [sum(abs.(resonant_flux[n_rows, s])) for s in 1:msing] bar!(p, 1:msing, rf_n; label="n=$nn", alpha=0.7) @@ -187,9 +188,9 @@ function plot_pe_delta_prime(h5path; save_path=nothing) _has_pe_data(h5path, key) || return plot(; title="No PE Δ' data — run with perturbed equilibrium enabled", legend=false) - delta_prime, q_sing, msing, mn_index = h5open(h5path, "r") do fid + delta_prime, q_sing, msing, pe_n = h5open(h5path, "r") do fid read(fid[key]), read(fid["singular/q"]), read(fid["singular/msing"]), - read(fid["info/mn_index"]) + read(fid["perturbed_equilibrium/forcing_modes/n"]) end labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] @@ -197,9 +198,9 @@ function plot_pe_delta_prime(h5path; save_path=nothing) p = plot(; xlabel="rational surface", ylabel="|Δ'|", title="Tearing stability Δ' (PE)", legend=:outertopright) - n_vals = unique(mn_index[:, 2]) + n_vals = unique(pe_n) for nn in n_vals - n_rows = findall(j -> mn_index[j, 2] == nn, 1:size(mn_index, 1)) + n_rows = findall(==(nn), pe_n) dp_n = [maximum(abs.(delta_prime[n_rows, s])) for s in 1:msing] bar!(p, 1:msing, dp_n; label="n=$nn", alpha=0.7) end @@ -253,9 +254,9 @@ function _plot_resonant_current(h5path) _has_pe_data(h5path, key) || return plot(; title="No resonant current data", legend=false) - resonant_current, q_sing, msing, mn_index = h5open(h5path, "r") do fid + resonant_current, q_sing, msing, pe_n = h5open(h5path, "r") do fid read(fid[key]), read(fid["singular/q"]), read(fid["singular/msing"]), - read(fid["info/mn_index"]) + read(fid["perturbed_equilibrium/forcing_modes/n"]) end labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] @@ -263,9 +264,9 @@ function _plot_resonant_current(h5path) p = plot(; xlabel="rational surface", ylabel="|I_res|", title="Resonant current |I_res| per surface", legend=:outertopright) - n_vals = unique(mn_index[:, 2]) + n_vals = unique(pe_n) for nn in n_vals - n_rows = findall(j -> mn_index[j, 2] == nn, 1:size(mn_index, 1)) + n_rows = findall(==(nn), pe_n) rc_n = [sum(abs.(resonant_current[n_rows, s])) for s in 1:msing] bar!(p, 1:msing, rc_n; label="n=$nn", alpha=0.7) end @@ -315,19 +316,19 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) _has_pe_data(h5path, base * rkey) || return plot(; title="No response data — run with perturbed equilibrium enabled", legend=false) - data_r, data_i, xs, mlow, mhigh, msing, psi_sing = h5open(h5path, "r") do fid + data_r, data_i, psi_response, mlow, mhigh, msing, psi_sing = h5open(h5path, "r") do fid read(fid[base * rkey]), read(fid[base * ikey]), - read(fid["splines/profiles/xs"]), + read(fid["integration/psi"]), read(fid["info/mlow"]), read(fid["info/mhigh"]), read(fid["singular/msing"]), read(fid["singular/psi"]) end - data = complex.(data_r, data_i) # shape: (numpert_total, npsi) + data = complex.(data_r, data_i) # shape: (npsi, numpert_total) mpert = mhigh - mlow + 1 # Use only the first n's modes for a clean spectrogram (single-n assumption for display) m_vals = mlow:mhigh - data_mn = data[1:mpert, :] # first mpert rows correspond to first n + data_mn = data[:, 1:mpert] # (npsi, mpert): first mpert columns = first n's modes # Top panel: line plot per mode p1 = plot(; @@ -338,12 +339,12 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) ) cmap = cgrad(:roma, mpert; categorical=true) for (i, m) in enumerate(m_vals) - plot!(p1, xs, abs.(data_mn[i, :]); label="m=$m", color=cmap[i], linewidth=1.5) + plot!(p1, psi_response, abs.(data_mn[:, i]); label="m=$m", color=cmap[i], linewidth=1.5) end - # Bottom panel: heatmap + # Bottom panel: heatmap — z must be (n_m_vals, n_psi) = (mpert, npsi) p2 = heatmap( - xs, collect(m_vals), abs.(data_mn); + psi_response, collect(m_vals), abs.(data_mn'); xlabel="ψ_N", ylabel="m", title="", @@ -365,8 +366,7 @@ end Three-panel composite summary of perturbed equilibrium results: - - Top-left: Island half-widths and Chirikov parameter overlay (`plot_island_widths` + - `plot_chirikov_parameter` on shared axes) + - Top-left: Island half-widths (`plot_island_widths`) - Top-right: Energy breakdown — plasma, vacuum, and total energies - Bottom: ξ_ψ mode spectrogram (`plot_mode_spectrogram`) From faf213aa807e3f32499332f961c850ca423200a3 Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Mon, 23 Mar 2026 12:31:34 -0400 Subject: [PATCH 03/12] ANALYSIS - IMPROVEMENT - Spectrogram axis convention, legend trimming, and spacing fixes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - plot_mode_spectrogram: flip heatmap to (m, ψ_N) with psi on vertical axis (conventional spectrogram orientation); rational surface overlays use hline! - plot_mode_spectrogram: trim line plot legend to m ∈ [0, nhigh·q95) to avoid unwieldy legend from negative/high-m modes outside the resonant range - All plots: add left_margin, bottom_margin, right_margin (Plots.mm) throughout to prevent axis labels, rotated xtick text, and colorbar titles from clipping - plot_equilibrium_summary: increase size to (1300, 450), add top_margin for plot_title clearance; reduce plot_titlefontsize to avoid overlap - Summary plots: adjust sizes to give subpanels more breathing room Co-Authored-By: Claude Sonnet 4.6 --- src/Analysis/Equilibrium.jl | 16 +++++-- src/Analysis/ForceFreeStates.jl | 46 +++++++++++++------ src/Analysis/PerturbedEquilibrium.jl | 68 ++++++++++++++++++---------- 3 files changed, 89 insertions(+), 41 deletions(-) diff --git a/src/Analysis/Equilibrium.jl b/src/Analysis/Equilibrium.jl index 12380dae8..ecc8ce641 100644 --- a/src/Analysis/Equilibrium.jl +++ b/src/Analysis/Equilibrium.jl @@ -97,7 +97,10 @@ function plot_qprofile(h5path; show_singular=true, save_path=nothing) xlabel="ψ_N", ylabel="q", title="Safety factor q(ψ)", - legend=false + legend=false, + left_margin=5Plots.mm, + bottom_margin=5Plots.mm, + right_margin=8Plots.mm # annotations at x=1.0 need room on the right ) hline!(p, [q0, q95]; linestyle=:dot, color=:gray, label=nothing) annotate!(p, 1.0, q0, text(" q0=$(round(q0, digits=2))", 7, :left, :gray)) @@ -147,7 +150,9 @@ function plot_pressure_profile(h5path; save_path=nothing) xlabel="ψ_N", ylabel="μ₀p", title="Pressure profile (βₜ=$(round(betat, digits=3)), βₚ=$(round(betap1, digits=3)))", - legend=false + legend=false, + left_margin=5Plots.mm, + bottom_margin=5Plots.mm ) for s in 1:msing vline!(p, [psi_sing[s]]; linestyle=:dash, color=:red, label=nothing) @@ -184,7 +189,9 @@ function plot_f_profile(h5path; save_path=nothing) xlabel="ψ_N", ylabel="2πF", title="Toroidal field function 2πF(ψ)", - legend=false + legend=false, + left_margin=5Plots.mm, + bottom_margin=5Plots.mm ) isnothing(save_path) || savefig(p, save_path) @@ -225,7 +232,8 @@ function plot_equilibrium_summary(h5path; save_path=nothing) p_pres = plot_pressure_profile(h5path) p_f = plot_f_profile(h5path) - p = plot(p_q, p_pres, p_f; layout=(1, 3), size=(1200, 400), plot_title=title_str) + p = plot(p_q, p_pres, p_f; layout=(1, 3), size=(1300, 450), + plot_title=title_str, plot_titlefontsize=10, top_margin=8Plots.mm) isnothing(save_path) || savefig(p, save_path) return p diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index 89e53758b..bdde5eb52 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -39,7 +39,9 @@ function plot_mode_displacement(h5path; modes=1:5, save_path=nothing) p = plot(; xlabel="ψ_N", ylabel="Im(ξ_ψ)", - title="Least Stable Eigenmode ξ_ψ" + title="Least Stable Eigenmode ξ_ψ", + left_margin=5Plots.mm, + bottom_margin=5Plots.mm ) for m in modes mlow <= m <= mhigh || continue @@ -96,14 +98,19 @@ function plot_eigenmode_summary(h5path; save_path=nothing) xlabel="m", ylabel="|Eigenvector|", title="Mode 1, |λ₁| = $(round(abs(et[1]), digits=3))", - legend=false + legend=false, + left_margin=5Plots.mm, + bottom_margin=5Plots.mm ) p2 = heatmap( m_vals, 1:nmodes, abs.(wt'); xlabel="m", ylabel="mode index", - colorbar_title="|Wₜ|" + colorbar_title="|Wₜ|", + left_margin=5Plots.mm, + right_margin=10Plots.mm, + bottom_margin=5Plots.mm ) p3 = scatter( @@ -111,11 +118,13 @@ function plot_eigenmode_summary(h5path; save_path=nothing) xlabel="|Eigenvalue|", ylabel="mode index", xscale=:log10, - legend=false + legend=false, + left_margin=5Plots.mm, + bottom_margin=5Plots.mm ) l = @layout [a{0.25h}; b c{0.25w}] - p = plot(p1, p2, p3; layout=l, size=(900, 700)) + p = plot(p1, p2, p3; layout=l, size=(950, 750)) isnothing(save_path) || savefig(p, save_path) return p @@ -149,7 +158,9 @@ function plot_stability_criterion(h5path; save_path=nothing) xlabel="ψ_N", ylabel="crit", title="Stability criterion (smallest eigenvalue of W⁻¹) vs ψ_N", - legend=false + legend=false, + left_margin=5Plots.mm, + bottom_margin=5Plots.mm ) isnothing(save_path) || savefig(p, save_path) @@ -201,7 +212,10 @@ function plot_energy_eigenvectors(h5path; matrix_type=:total, save_path=nothing) xlabel="m", ylabel="mode index", title="Energy eigenvectors |Wₜ| (total)", - colorbar_title="|Wₜ|" + colorbar_title="|Wₜ|", + left_margin=5Plots.mm, + right_margin=10Plots.mm, + bottom_margin=5Plots.mm ) isnothing(save_path) || savefig(p, save_path) @@ -248,7 +262,9 @@ function plot_eigenvalue_spectrum(h5path; matrix_type=:total, save_path=nothing) title="Eigenvalue spectrum ($matrix_type)", legend=false, color=colors, - markerstrokewidth=0 + markerstrokewidth=0, + left_margin=5Plots.mm, + bottom_margin=5Plots.mm ) vline!(p, [0]; linestyle=:dash, color=:black, label=nothing) @@ -319,7 +335,9 @@ function plot_delta_prime(h5path; save_path=nothing) title="Tearing stability Δ' (FFS asymptotic coefficients)", legend=false, color=colors, - xrotation=30 + xrotation=30, + bottom_margin=10Plots.mm, + left_margin=5Plots.mm ) isnothing(save_path) || savefig(p, save_path) @@ -362,10 +380,10 @@ function plot_dcon_summary(h5path; save_path=nothing) if has_vac p_evec = plot_energy_eigenvectors(h5path; matrix_type=:total) p_eval = plot_eigenvalue_spectrum(h5path; matrix_type=:total) - p = plot(p_evec, p_crit, p_eval, p_dp; layout=(2, 2), size=(1000, 800)) + p = plot(p_evec, p_crit, p_eval, p_dp; layout=(2, 2), size=(1100, 900)) else title!(p_crit, "Stability criterion (no vacuum data — rerun with vac_flag = true)") - p = plot(p_crit, p_dp; layout=(1, 2), size=(1000, 400)) + p = plot(p_crit, p_dp; layout=(1, 2), size=(1100, 500)) end isnothing(save_path) || savefig(p, save_path) @@ -404,7 +422,9 @@ function plot_singular_surfaces(h5path; save_path=nothing) xlabel="ψ_N", ylabel="q", title="Safety factor and rational surfaces", - legend=false + legend=false, + left_margin=5Plots.mm, + bottom_margin=5Plots.mm ) hline!(p_q, [q0, q95]; linestyle=:dot, color=:gray, label=nothing) for s in 1:msing @@ -415,7 +435,7 @@ function plot_singular_surfaces(h5path; save_path=nothing) p_dp = plot_delta_prime(h5path) - p = plot(p_q, p_dp; layout=(1, 2), size=(1000, 400)) + p = plot(p_q, p_dp; layout=(1, 2), size=(1100, 500)) isnothing(save_path) || savefig(p, save_path) return p diff --git a/src/Analysis/PerturbedEquilibrium.jl b/src/Analysis/PerturbedEquilibrium.jl index f77a06d5e..5f5b24468 100644 --- a/src/Analysis/PerturbedEquilibrium.jl +++ b/src/Analysis/PerturbedEquilibrium.jl @@ -52,7 +52,8 @@ function plot_resonant_flux(h5path; save_path=nothing) xticks_arg = (1:msing, labels) p = plot(; xlabel="rational surface", ylabel="|Φ_res|", - title="Resonant flux |Φ_res| per surface", legend=:outertopright) + title="Resonant flux |Φ_res| per surface", legend=:outertopright, + left_margin=5Plots.mm, bottom_margin=10Plots.mm) n_vals = unique(pe_n) for nn in n_vals @@ -105,7 +106,9 @@ function plot_island_widths(h5path; save_path=nothing) title="Island half-widths", legend=false, color=:steelblue, - xrotation=30 + xrotation=30, + left_margin=5Plots.mm, + bottom_margin=10Plots.mm ) isnothing(save_path) || savefig(p, save_path) @@ -152,7 +155,9 @@ function plot_chirikov_parameter(h5path; save_path=nothing) title="Chirikov overlap parameter", legend=false, color=colors, - xrotation=30 + xrotation=30, + left_margin=5Plots.mm, + bottom_margin=10Plots.mm ) hline!(p, [1.0]; linestyle=:dash, color=:black, label=nothing) @@ -196,7 +201,8 @@ function plot_pe_delta_prime(h5path; save_path=nothing) labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] p = plot(; xlabel="rational surface", ylabel="|Δ'|", - title="Tearing stability Δ' (PE)", legend=:outertopright) + title="Tearing stability Δ' (PE)", legend=:outertopright, + left_margin=5Plots.mm, bottom_margin=10Plots.mm) n_vals = unique(pe_n) for nn in n_vals @@ -242,7 +248,7 @@ function plot_resonant_field(h5path; save_path=nothing) p4 = plot_island_widths(h5path) p5 = plot_chirikov_parameter(h5path) - p = plot(p1, p2, p3, p4, p5; layout=(5, 1), size=(800, 1400)) + p = plot(p1, p2, p3, p4, p5; layout=(5, 1), size=(900, 1600)) isnothing(save_path) || savefig(p, save_path) return p @@ -262,7 +268,8 @@ function _plot_resonant_current(h5path) labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] p = plot(; xlabel="rational surface", ylabel="|I_res|", - title="Resonant current |I_res| per surface", legend=:outertopright) + title="Resonant current |I_res| per surface", legend=:outertopright, + left_margin=5Plots.mm, bottom_margin=10Plots.mm) n_vals = unique(pe_n) for nn in n_vals @@ -280,9 +287,10 @@ end Two-panel spectrogram of a perturbed equilibrium response field component: - - Top: `|component|` vs ψ_N, one curve per poloidal mode m (colored by m) - - Bottom: Heatmap of `|component|` in (ψ_N × m) space, with white dashed lines at - rational surface locations + - Top: `|component|` vs ψ_N, one curve per poloidal mode m. Only resonant modes + (m ∈ [0, nhigh·q95)) are labeled to keep the legend readable. + - Bottom: Heatmap of `|component|` in (m, ψ_N) space (psi on vertical axis), with + white dashed lines at rational surface locations. Inspired by `plot_spectrograms.py` from OMFIT GPEC. @@ -316,10 +324,11 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) _has_pe_data(h5path, base * rkey) || return plot(; title="No response data — run with perturbed equilibrium enabled", legend=false) - data_r, data_i, psi_response, mlow, mhigh, msing, psi_sing = h5open(h5path, "r") do fid + data_r, data_i, psi_response, mlow, mhigh, nhigh, q95, msing, psi_sing = h5open(h5path, "r") do fid read(fid[base * rkey]), read(fid[base * ikey]), read(fid["integration/psi"]), - read(fid["info/mlow"]), read(fid["info/mhigh"]), + read(fid["info/mlow"]), read(fid["info/mhigh"]), read(fid["info/nhigh"]), + read(fid["equil/q95"]), read(fid["singular/msing"]), read(fid["singular/psi"]) end @@ -330,32 +339,41 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) m_vals = mlow:mhigh data_mn = data[:, 1:mpert] # (npsi, mpert): first mpert columns = first n's modes - # Top panel: line plot per mode + # Top panel: line plot per mode — only label resonant range m ∈ [0, nhigh·q95) + m_max_legend = nhigh * q95 p1 = plot(; xlabel="ψ_N", ylabel="|$(component)|", title="Mode spectrogram: $(component)", - legend=:outertopright + legend=:outertopright, + left_margin=5Plots.mm, + bottom_margin=5Plots.mm ) cmap = cgrad(:roma, mpert; categorical=true) for (i, m) in enumerate(m_vals) - plot!(p1, psi_response, abs.(data_mn[:, i]); label="m=$m", color=cmap[i], linewidth=1.5) + show_label = 0 <= m < m_max_legend + plot!(p1, psi_response, abs.(data_mn[:, i]); + label=(show_label ? "m=$m" : nothing), color=cmap[i], linewidth=1.5) end - # Bottom panel: heatmap — z must be (n_m_vals, n_psi) = (mpert, npsi) + # Bottom panel: heatmap in (m, ψ_N) space — psi on vertical axis, m on horizontal + # z must be (n_psi, n_m) = (npsi, mpert) so that rows = psi, cols = m p2 = heatmap( - psi_response, collect(m_vals), abs.(data_mn'); - xlabel="ψ_N", - ylabel="m", + collect(m_vals), psi_response, abs.(data_mn); + xlabel="m", + ylabel="ψ_N", title="", - colorbar_title="|$(component)|" + colorbar_title="|$(component)|", + left_margin=5Plots.mm, + right_margin=10Plots.mm, + bottom_margin=5Plots.mm ) # Overlay rational surface locations as white dashed lines for s in 1:msing - vline!(p2, [psi_sing[s]]; linestyle=:dash, color=:white, linewidth=1.5, label=nothing) + hline!(p2, [psi_sing[s]]; linestyle=:dash, color=:white, linewidth=1.5, label=nothing) end - p = plot(p1, p2; layout=(2, 1), size=(900, 700)) + p = plot(p1, p2; layout=(2, 1), size=(950, 750)) isnothing(save_path) || savefig(p, save_path) return p @@ -388,7 +406,7 @@ function plot_perturbed_equilibrium_summary(h5path; save_path=nothing) p_spectro = plot_mode_spectrogram(h5path; component=:xi_psi) l = @layout [grid(1, 2){0.35h}; b] - p = plot(p_islands, p_energies, p_spectro; layout=l, size=(1000, 1000)) + p = plot(p_islands, p_energies, p_spectro; layout=l, size=(1100, 1100)) isnothing(save_path) || savefig(p, save_path) return p @@ -406,7 +424,7 @@ function _plot_energies(h5path) read(fid["perturbed_equilibrium/energies/total_energy"]) end - vals = [real(ep), real(ev), real(et)] + vals = [real(ep), real(ev), real(et)] # real() in case energies are stored as complex names = ["Plasma", "Vacuum", "Total"] colors = [:steelblue, :darkorange, :green] @@ -415,7 +433,9 @@ function _plot_energies(h5path) ylabel="Energy", title="Energy breakdown", legend=false, - color=colors + color=colors, + left_margin=5Plots.mm, + bottom_margin=5Plots.mm ) hline!(p, [0]; linestyle=:dash, color=:black, label=nothing) From 910df8f17d63fcba024f2757810421e8e9ef83b5 Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Mon, 23 Mar 2026 12:35:12 -0400 Subject: [PATCH 04/12] ANALYSIS - IMPROVEMENT - Fix eigenmode summary layout and eigenvalue coloring - plot_eigenmode_summary: align top eigenvector bar with heatmap x-axis using @layout [a{0.25h} _{0.25w}; b c{0.25w}] (blank cell preserves column widths) - plot_eigenmode_summary: color scatter by sign (red=negative, blue=positive) and add markerstrokewidth=0 to remove black rings from dots - plot_eigenvalue_spectrum: apply same color convention (red=negative/unstable, blue=positive/stable) instead of previous red/green Co-Authored-By: Claude Sonnet 4.6 --- src/Analysis/ForceFreeStates.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index bdde5eb52..040599cc0 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -113,17 +113,21 @@ function plot_eigenmode_summary(h5path; save_path=nothing) bottom_margin=5Plots.mm ) + colors = [real(e) < 0 ? :red : :blue for e in et] p3 = scatter( abs.(et), 1:nmodes; xlabel="|Eigenvalue|", ylabel="mode index", xscale=:log10, legend=false, + color=colors, + markerstrokewidth=0, left_margin=5Plots.mm, bottom_margin=5Plots.mm ) - l = @layout [a{0.25h}; b c{0.25w}] + # Top panel (p1) shares m-axis with heatmap (p2); blank cell keeps widths aligned + l = @layout [a{0.25h} _{0.25w}; b c{0.25w}] p = plot(p1, p2, p3; layout=l, size=(950, 750)) isnothing(save_path) || savefig(p, save_path) @@ -253,7 +257,7 @@ function plot_eigenvalue_spectrum(h5path; matrix_type=:total, save_path=nothing) nmodes = length(et) ev_real = real.(et) - colors = [v > 0 ? :red : :green for v in ev_real] + colors = [v < 0 ? :red : :blue for v in ev_real] # red = negative (unstable), blue = positive (stable) p = scatter( ev_real, 1:nmodes; From 66b199bcaf4319aee6eba3d6505d45e39ec91a0a Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Wed, 25 Mar 2026 17:35:40 -0400 Subject: [PATCH 05/12] ANALYSIS - IMPROVEMENT - Polish plots: layout, labels, scatter/line conversions, new equilibrium panels MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Eigenmode summary (ffs_eigenmode_summary): removed from make_plots (redundant with ffs_summary) - ffs_dcon_summary → ffs_summary; eigenvalue plot replaced with mode displacement (|ξ_ψ| vs ψ_N) - All bar charts (Δ', resonant flux/current, island widths, Chirikov) converted to scatter plots on ψ_N axis with integer-q annotations - plot_delta_prime: shows Re(Δ') with sign coloring and zero reference line - Equilibrium summary: tall RZ panel (left column), profiles stacked right; x-axes share [0,1] range; upper panels suppress redundant x-axis labels - plot_flux_surfaces_h5: rational surfaces drawn as red closed contours with q labels - plot_f_profile: rational surface vlines added (consistent with q and P plots) - GSE panels: log-scale y-axis; integrated error overplotted on theta-slice panel; single legend entry - All ψ_N labels use LaTeXStrings (L"\psi_N"); |ξ_ψ| and |D_c| similarly subscripted - Profile plot titles removed (y-axis labels carry the information) - Stability criterion: ylabel/title updated to |D_c|, zero reference line added - PlotlyJS and LaTeXStrings added to Project.toml - CLAUDE.md: note added to never remove packages from Project.toml Co-Authored-By: Claude Sonnet 4.6 --- CLAUDE.md | 1 + Project.toml | 4 + examples/DIIID-like_ideal_example/gpec.toml | 108 ++++---- .../DIIID-like_ideal_example/make_plots.jl | 10 +- src/Analysis/Equilibrium.jl | 240 +++++++++++++++++- src/Analysis/ForceFreeStates.jl | 92 +++---- src/Analysis/PerturbedEquilibrium.jl | 139 +++++----- 7 files changed, 420 insertions(+), 174 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 7d7c8c9d3..d88b49ab3 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -475,6 +475,7 @@ This format is used for compiling release notes, so tags should be human-readabl ### General - **Julia version**: 1.11 is the target version +- **Never remove packages from Project.toml** - If a package fails to load or resolve, run `Pkg.add(...)` or `Pkg.instantiate()` to fix the local environment. Do NOT remove the package from `Project.toml`. The developer works across multiple branches and machines, so environment drift is expected — the right fix is always to update the environment to satisfy the toml, not to trim the toml to match the current environment state. - **Indexing**: The codebase uses 0-based indexing in many places to match Fortran conventions, then converts to 1-based Julia indexing - **No step numbering in code comments** - Avoid annotations like "Step 1: do this" followed by "Step 2: do that". These get out of sync as code changes. Just describe the action without numbering. - **Keep code comments concise** - A comment should be one line where possible. Do not write multi-line block comments explaining the current session's investigation, what was tried, what was wrong before, or why a specific file/path behaves differently. State what the code does and why at a general level. Example of too much detail: a 6-line block explaining that efit_by_inversion uses psilow>0 while CHEASE starts at 0, that the old code was removed, and that spline spikes result. Preferred: `# Replicate Fortran inverse.f: overwrite deta at axis (r²=0) by extrapolating from innermost surfaces.` diff --git a/Project.toml b/Project.toml index 65c64c3ae..9e480fa04 100644 --- a/Project.toml +++ b/Project.toml @@ -15,9 +15,11 @@ FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FastInterpolations = "9ea80cae-fc13-4c00-8066-6eaedb12f34b" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" @@ -39,9 +41,11 @@ FastGaussQuadrature = "1.1.0" FastInterpolations = "0.2.13" HDF5 = "0.17.2" JLD2 = "0.6.3" +LaTeXStrings = "1.4.0" LinearAlgebra = "1" OrdinaryDiffEq = "6.102.0" Pkg = "1" +PlotlyJS = "0.18.17" Plots = "1.40.15" Printf = "1" Roots = "2.2.13" diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 7e8cf0216..25c660829 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -1,68 +1,68 @@ [Equilibrium] -eq_type = "efit" # Type of the input 2D equilibrium file eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file -jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc) -power_bp = 0 # Poloidal field power exponent for Jacobian -power_b = 0 # Toroidal field power exponent for Jacobian -power_r = 0 # Major radius power exponent for Jacobian -grid_type = "log_asymptotic" # Radial grid packing type -psilow = 1e-4 # Lower limit of normalized flux coordinate -psihigh = 0.993 # Upper limit of normalized flux coordinate -mpsi = 0 # Number of radial grid points (0 = auto-compute from psi_accuracy) -psi_accuracy = 0.001 # Target absolute error in q for auto-mpsi -mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) -etol = 1e-7 # Error tolerance for equilibrium solver -force_termination = false # Terminate after equilibrium setup (skip stability calculations) +eq_type = "efit" # Type of the input 2D equilibrium file +jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc) +power_bp = 0 # Poloidal field power exponent for Jacobian +power_b = 0 # Toroidal field power exponent for Jacobian +power_r = 0 # Major radius power exponent for Jacobian +grid_type = "log_asymptotic" # Radial grid packing type +psilow = 1e-4 # Lower limit of normalized flux coordinate +psihigh = 0.993 # Upper limit of normalized flux coordinate +mpsi = 0 # Number of radial grid points (0 = auto-compute from psi_accuracy) +psi_accuracy = 0.001 # Target absolute error in q for auto-mpsi +mtheta = 256 # Number of poloidal grid points +newq0 = 0 # Override for on-axis safety factor (0 = use input value) +etol = 1e-7 # Error tolerance for equilibrium solver +force_termination = false # Terminate after equilibrium setup (skip stability calculations) [Wall] -shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) -equal_arc_wall = true # Equal arc length distribution of nodes on wall +shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) +a = 0.2415 # Distance from plasma (conformal) or shape parameter +aw = 0.05 # Half-thickness parameter for Dee-shaped walls +bw = 1.5 # Elongation parameter for wall shapes +cw = 0 # Offset of wall center from major radius +dw = 0.5 # Triangularity parameter for wall shapes +tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForceFreeStates] -bal_flag = false # Ideal MHD ballooning criterion for short wavelengths -mat_flag = true # Construct coefficient matrices for diagnostic purposes -ode_flag = true # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC) -vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes -mer_flag = true # Evaluate the Mercier criterian +bal_flag = false # Ideal MHD ballooning criterion for short wavelengths +mat_flag = true # Construct coefficient matrices for diagnostic purposes +ode_flag = true # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC) +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes +mer_flag = true # Evaluate the Mercier criterian -set_psilim_via_dmlim = true # Safety factor (q) limit determined as q_ir+dmlim... -dmlim = 0.2 # See set_psilim_via_dmlim -psiedge = 1.00 # If less then psilim, calculates dW(psi)... -qlow = 1.02 # Integration initiated at q determined by min(q0, qlow)... -qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh)... -sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) +set_psilim_via_dmlim = true # Safety factor (q) limit determined as q_ir+dmlim... +dmlim = 0.2 # See set_psilim_via_dmlim +psiedge = 1.00 # If less then psilim, calculates dW(psi)... +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow)... +qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh)... +sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) -nn_low = 1 # Smallest toroidal mode number to include -nn_high = 1 # Largest toroidal mode number to include -delta_mlow = 8 # Expands lower bound of Fourier harmonics -delta_mhigh = 8 # Expands upper bound of Fourier harmonics -delta_mband = 0 # Integration keeps only this wide a band... -mthvac = 512 # Number of points used in splines over poloidal angle at plasma-vacuum interface. -thmax0 = 1 # Linear multiplier on the automatic choice of theta integration bounds +nn_low = 1 # Smallest toroidal mode number to include +nn_high = 1 # Largest toroidal mode number to include +delta_mlow = 8 # Expands lower bound of Fourier harmonics +delta_mhigh = 8 # Expands upper bound of Fourier harmonics +delta_mband = 0 # Integration keeps only this wide a band... +mthvac = 512 # Number of points used in splines over poloidal angle at plasma-vacuum interface. +thmax0 = 1 # Linear multiplier on the automatic choice of theta integration bounds -kin_flag = false # Kinetic EL equation (default: false) -con_flag = false # Continue integration through layers (default: false) -kinfac1 = 1.0 # Scale factor for energy contribution (default: 1.0) -kinfac2 = 1.0 # Scale factor for torque contribution (default: 1.0) -kingridtype = 0 # Regular grid method (default: 0) -passing_flag = true # Includes passing particle effects (default: false) -ktanh_flag = true # Ignore kinetic effects in the core smoothly (default: false) -ktc = 0.1 # Parameter for ktanh_flag (default: 0.1) -ktw = 50.0 # Parameter for ktanh_flag (default: 50.0) -ion_flag = true # Include ion dW_k when kin_flag is true -electron_flag = false # Include electron dW_k when kin_flag is true +kin_flag = false # Kinetic EL equation (default: false) +con_flag = false # Continue integration through layers (default: false) +kinfac1 = 1.0 # Scale factor for energy contribution (default: 1.0) +kinfac2 = 1.0 # Scale factor for torque contribution (default: 1.0) +kingridtype = 0 # Regular grid method (default: 0) +passing_flag = true # Includes passing particle effects (default: false) +ktanh_flag = true # Ignore kinetic effects in the core smoothly (default: false) +ktc = 0.1 # Parameter for ktanh_flag (default: 0.1) +ktw = 50.0 # Parameter for ktanh_flag (default: 50.0) +ion_flag = true # Include ion dW_k when kin_flag is true +electron_flag = false # Include electron dW_k when kin_flag is true eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations -save_interval = 3 # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. -singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e4 # Maximum fraction of solutions allowed before re-normalized +save_interval = 1 # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. +singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced +ucrit = 1e4 # Maximum fraction of solutions allowed before re-normalized [ForcingTerms] forcing_data_file = "forcing.dat" # Path to forcing data file (n, m, complex amplitude) diff --git a/examples/DIIID-like_ideal_example/make_plots.jl b/examples/DIIID-like_ideal_example/make_plots.jl index 63c6f9f66..68e728e0a 100644 --- a/examples/DIIID-like_ideal_example/make_plots.jl +++ b/examples/DIIID-like_ideal_example/make_plots.jl @@ -11,6 +11,9 @@ mkpath(outdir) p = Analysis.Equilibrium.plot_equilibrium_summary(h5path) savefig(p, joinpath(outdir, "equil_summary.png")) +p = Analysis.Equilibrium.plot_flux_surfaces_h5(h5path) +savefig(p, joinpath(outdir, "equil_flux_surfaces.png")) + p = Analysis.Equilibrium.plot_qprofile(h5path) savefig(p, joinpath(outdir, "equil_qprofile.png")) @@ -24,9 +27,6 @@ savefig(p, joinpath(outdir, "equil_F.png")) p = Analysis.ForceFreeStates.plot_mode_displacement(h5path; modes=1:5) savefig(p, joinpath(outdir, "ffs_mode_displacement.png")) -p = Analysis.ForceFreeStates.plot_eigenmode_summary(h5path) -savefig(p, joinpath(outdir, "ffs_eigenmode_summary.png")) - p = Analysis.ForceFreeStates.plot_stability_criterion(h5path) savefig(p, joinpath(outdir, "ffs_stability_criterion.png")) @@ -39,8 +39,8 @@ savefig(p, joinpath(outdir, "ffs_eigenvalue_spectrum.png")) p = Analysis.ForceFreeStates.plot_delta_prime(h5path) savefig(p, joinpath(outdir, "ffs_delta_prime.png")) -p = Analysis.ForceFreeStates.plot_dcon_summary(h5path) -savefig(p, joinpath(outdir, "ffs_dcon_summary.png")) +p = Analysis.ForceFreeStates.plot_ffs_summary(h5path) +savefig(p, joinpath(outdir, "ffs_summary.png")) p = Analysis.ForceFreeStates.plot_singular_surfaces(h5path) savefig(p, joinpath(outdir, "ffs_singular_surfaces.png")) diff --git a/src/Analysis/Equilibrium.jl b/src/Analysis/Equilibrium.jl index ecc8ce641..46cf3b448 100644 --- a/src/Analysis/Equilibrium.jl +++ b/src/Analysis/Equilibrium.jl @@ -6,6 +6,7 @@ Post-processing and visualization functions for GPEC equilibrium objects and HDF module Equilibrium using HDF5 +using LaTeXStrings using Plots """ @@ -94,10 +95,11 @@ function plot_qprofile(h5path; show_singular=true, save_path=nothing) p = plot( xs, q; - xlabel="ψ_N", + xlabel=L"\psi_N", ylabel="q", - title="Safety factor q(ψ)", + title="", legend=false, + xlims=(0, 1), left_margin=5Plots.mm, bottom_margin=5Plots.mm, right_margin=8Plots.mm # annotations at x=1.0 need room on the right @@ -139,18 +141,18 @@ 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, betat, betap1, msing, psi_sing = h5open(h5path, "r") do fid + xs, mu0p, msing, psi_sing = h5open(h5path, "r") do fid read(fid["splines/profiles/xs"]), read(fid["splines/profiles/mu0p"]), - read(fid["equil/betat"]), read(fid["equil/betap1"]), read(fid["singular/msing"]), read(fid["singular/psi"]) end p = plot( xs, mu0p; - xlabel="ψ_N", + xlabel=L"\psi_N", ylabel="μ₀p", - title="Pressure profile (βₜ=$(round(betat, digits=3)), βₚ=$(round(betap1, digits=3)))", + title="", legend=false, + xlims=(0, 1), left_margin=5Plots.mm, bottom_margin=5Plots.mm ) @@ -180,15 +182,206 @@ 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 = h5open(h5path, "r") do fid - read(fid["splines/profiles/xs"]), read(fid["splines/profiles/2piF"]) + xs, twopif, msing, psi_sing = h5open(h5path, "r") do fid + read(fid["splines/profiles/xs"]), read(fid["splines/profiles/2piF"]), + read(fid["singular/msing"]), read(fid["singular/psi"]) end p = plot( xs, twopif; - xlabel="ψ_N", + xlabel=L"\psi_N", ylabel="2πF", - title="Toroidal field function 2πF(ψ)", + title="", + legend=false, + xlims=(0, 1), + left_margin=5Plots.mm, + bottom_margin=5Plots.mm + ) + for s in 1:msing + vline!(p, [psi_sing[s]]; linestyle=:dash, color=:red, label=nothing) + end + + isnothing(save_path) || savefig(p, save_path) + return p +end + +""" + plot_flux_surfaces_h5(h5path; 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. + +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 + +### Keyword arguments + + - `n_psi`: Number of constant-ψ contours to draw (default: 11) + - `n_theta`: Number of constant-θ/θ spokes to draw (default: 18) + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object. +""" +function plot_flux_surfaces_h5(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["splines/rzphi/rcoords"]), read(fid["splines/rzphi/offset"]), + read(fid["splines/rzphi/xs"]), read(fid["splines/rzphi/ys"]), + read(fid["equil/ro"]), read(fid["equil/zo"]), + read(fid["singular/msing"]), read(fid["singular/psi"]), read(fid["singular/q"]) + end + + 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 + + p = plot(; + title="Flux surfaces in (R, Z)", + xlabel="R [m]", + ylabel="Z [m]", + aspect_ratio=:equal, + legend=:outertopright, + left_margin=5Plots.mm, + bottom_margin=5Plots.mm + ) + + psi_indices = round.(Int, range(1, n_psi_grid; length=n_psi)) + theta_indices = round.(Int, range(1, n_theta_grid; length=n_theta)) + + for (i, ipsi) in enumerate(psi_indices) + label = i == 1 ? "Constant ψ" : "" + plot!(p, [R_grid[ipsi, :]; R_grid[ipsi, 1]], [Z_grid[ipsi, :]; Z_grid[ipsi, 1]]; + color=:steelblue, linewidth=1.5, label=label) + end + for (i, itheta) in enumerate(theta_indices) + label = i == 1 ? "Constant θ" : "" + plot!(p, R_grid[:, itheta], Z_grid[:, itheta]; + color=:tomato, linewidth=0.8, label=label) + end + # Draw rational surface flux contours in red, annotated at the outboard midplane + for s in 1:msing + 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))" + plot!(p, [R_grid[idx, :]; R_grid[idx, 1]], [Z_grid[idx, :]; Z_grid[idx, 1]]; + color=:red, linewidth=1.5, + label=s == 1 ? "Rational surface" : "") + annotate!(p, R_grid[idx, 1], Z_grid[idx, 1], text(" $q_label", 7, :left, :red)) + end + + isnothing(save_path) || savefig(p, save_path) + return p +end + +""" + plot_gse_by_theta(h5path; n_theta_lines=8, save_path=nothing) + +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). + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file (used to locate `gse.h5`) + +### Keyword arguments + + - `n_theta_lines`: Number of evenly-spaced θ slices to overlay (default: 8) + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object, or `nothing` if `gse.h5` is absent. +""" +function plot_gse_by_theta(h5path; n_theta_lines=8, save_path=nothing) + 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 + end + + _, ntheta, _ = size(gse_data) + theta_indices = round.(Int, range(1, ntheta; length=n_theta_lines)) + + p = plot(; + xlabel=L"\psi_N", + ylabel="GSE error", + title="", + yscale=:log10, + xlims=(0, 1), + left_margin=5Plots.mm, + bottom_margin=5Plots.mm + ) + for (i, itheta) in enumerate(theta_indices) + psi_col = gse_data[:, itheta, 2] + err_col = gse_data[:, itheta, 7] + plot!(p, psi_col, err_col; label=i == 1 ? "Const. θ" : "", alpha=0.6) + end + + # Overplot flux-surface-integrated error as thick reference line + gsei_path = joinpath(dirname(h5path), "gsei.h5") + if isfile(gsei_path) + xs_int, errori = h5open(gsei_path, "r") do fid + read(fid["xs"]), read(fid["errori"]) + end + plot!(p, xs_int, vec(errori); + color=:black, linewidth=2.5, linestyle=:solid, label="integrated") + end + + isnothing(save_path) || savefig(p, save_path) + return p +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). + +### Arguments + + - `h5path`: Path to a GPEC HDF5 output file (used to locate `gsei.h5`) + +### Keyword arguments + + - `save_path`: If provided, save the figure to this path (default: `nothing`) + +### Returns + +A `Plots.jl` plot object, or `nothing` if `gsei.h5` is absent. +""" +function plot_gse_integrated(h5path; save_path=nothing) + 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"]) + end + + p = plot( + xs, vec(errlogi); + xlabel=L"\psi_N", + ylabel="log₁₀(integrated GSE)", + title="Flux-surface-integrated Grad-Shafranov error", legend=false, left_margin=5Plots.mm, bottom_margin=5Plots.mm @@ -201,12 +394,16 @@ end """ plot_equilibrium_summary(h5path; save_path=nothing) -Three-panel summary of equilibrium profiles: +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: - q(ψ) safety factor with rational surface markers (`plot_qprofile`) - μ₀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. + ### Arguments - `h5path`: Path to a GPEC HDF5 output file @@ -228,12 +425,27 @@ function plot_equilibrium_summary(h5path; save_path=nothing) 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_h5(h5path) p_q = plot_qprofile(h5path; show_singular=true) p_pres = plot_pressure_profile(h5path) p_f = plot_f_profile(h5path) - - p = plot(p_q, p_pres, p_f; layout=(1, 3), size=(1300, 450), - plot_title=title_str, plot_titlefontsize=10, top_margin=8Plots.mm) + p_gse = plot_gse_by_theta(h5path) # includes integrated overplot; nothing if absent + + # 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. + hide_xaxis!(p) = plot!(p; xlabel="", xformatter=_->"", bottom_margin=1Plots.mm) + + if isnothing(p_gse) + hide_xaxis!(p_q); hide_xaxis!(p_pres) + l = @layout [a{0.38w} [b; c; d]] + p = plot(p_rz, p_q, p_pres, p_f; layout=l, size=(1300, 750), + plot_title=title_str, plot_titlefontsize=10, top_margin=8Plots.mm) + else + hide_xaxis!(p_q); hide_xaxis!(p_pres); hide_xaxis!(p_f) + l = @layout [a{0.38w} [b; c; d; e]] + p = plot(p_rz, p_q, p_pres, p_f, p_gse; layout=l, size=(1300, 1000), + plot_title=title_str, plot_titlefontsize=10, top_margin=8Plots.mm) + end isnothing(save_path) || savefig(p, save_path) return p diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index 040599cc0..08683595a 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -7,13 +7,14 @@ results stored in GPEC HDF5 output files. module ForceFreeStates using HDF5 +using LaTeXStrings using Plots """ plot_mode_displacement(h5path; modes=1:5, save_path=nothing) -Plot Im(ξ_ψ) vs ψ_N for the least stable eigenmode, showing one curve per requested -poloidal mode number m. +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 @@ -29,23 +30,27 @@ poloidal mode number m. A `Plots.jl` plot object. """ function plot_mode_displacement(h5path; modes=1:5, save_path=nothing) - mlow, xi_psi, psi = h5open(h5path, "r") do fid - read(fid["info/mlow"]), read(fid["integration/xi_psi"]), read(fid["integration/psi"]) + mlow, xi_psi, psi, et = h5open(h5path, "r") do fid + read(fid["info/mlow"]), read(fid["integration/xi_psi"]), + read(fid["integration/psi"]), read(fid["vacuum/et"]) 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, dW = $(round(real(dW), sigdigits=4))" p = plot(; - xlabel="ψ_N", - ylabel="Im(ξ_ψ)", - title="Least Stable Eigenmode ξ_ψ", + xlabel=L"\psi_N", + ylabel=L"|\xi_\psi|", + title=title_str, left_margin=5Plots.mm, bottom_margin=5Plots.mm ) for m in modes mlow <= m <= mhigh || continue - plot!(p, psi, imag.(xi_psi[m-mlow+1, 1, :]); label="m=$m") # DCON phase convention: ξ_ψ is purely imaginary on the real axis + plot!(p, psi, abs.(xi_psi[m-mlow+1, 1, :]); label="m=$m") end isnothing(save_path) || savefig(p, save_path) @@ -126,9 +131,11 @@ function plot_eigenmode_summary(h5path; save_path=nothing) bottom_margin=5Plots.mm ) - # Top panel (p1) shares m-axis with heatmap (p2); blank cell keeps widths aligned + # Top panel (p1) shares m-axis with heatmap (p2); blank cell keeps widths aligned. + # Plots.jl assigns plots sequentially to ALL layout slots including `_`, so pass an + # explicit empty plot() to fill the blank top-right cell. l = @layout [a{0.25h} _{0.25w}; b c{0.25w}] - p = plot(p1, p2, p3; layout=l, size=(950, 750)) + p = plot(p1, plot(), p2, p3; layout=l, size=(950, 750)) isnothing(save_path) || savefig(p, save_path) return p @@ -159,13 +166,14 @@ function plot_stability_criterion(h5path; save_path=nothing) p = plot( psi, crit; - xlabel="ψ_N", - ylabel="crit", - title="Stability criterion (smallest eigenvalue of W⁻¹) vs ψ_N", + xlabel=L"\psi_N", + ylabel=L"|D_c|", + title="Fixed-boundary stability", legend=false, left_margin=5Plots.mm, bottom_margin=5Plots.mm ) + hline!(p, [0.0]; linestyle=:dash, color=:black, label=nothing) isnothing(save_path) || savefig(p, save_path) return p @@ -260,9 +268,9 @@ function plot_eigenvalue_spectrum(h5path; matrix_type=:total, save_path=nothing) colors = [v < 0 ? :red : :blue for v in ev_real] # red = negative (unstable), blue = positive (stable) p = scatter( - ev_real, 1:nmodes; - xlabel="Re(eigenvalue)", - ylabel="mode index", + 1:nmodes, ev_real; + xlabel="mode index", + ylabel="Re(eigenvalue)", title="Eigenvalue spectrum ($matrix_type)", legend=false, color=colors, @@ -270,7 +278,7 @@ function plot_eigenvalue_spectrum(h5path; matrix_type=:total, save_path=nothing) left_margin=5Plots.mm, bottom_margin=5Plots.mm ) - vline!(p, [0]; linestyle=:dash, color=:black, label=nothing) + hline!(p, [0]; linestyle=:dash, color=:black, label=nothing) isnothing(save_path) || savefig(p, save_path) return p @@ -279,9 +287,9 @@ end """ plot_delta_prime(h5path; save_path=nothing) -Bar chart of `|Δ'|` per singular surface, computed from the stored asymptotic coefficients -`ca_left` and `ca_right`. Bars are colored red (tearing unstable, Re(Δ') > 0) or blue -(tearing stable). +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. Δ' 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`. @@ -299,8 +307,8 @@ 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, q_sing, ca_l, ca_r, psio, mn_index = h5open(h5path, "r") do fid - read(fid["singular/msing"]), read(fid["singular/q"]), + msing, psi_sing, q_sing, ca_l, ca_r, psio, mn_index = h5open(h5path, "r") do fid + read(fid["singular/msing"]), read(fid["singular/psi"]), read(fid["singular/q"]), read(fid["singular/ca_left"]), read(fid["singular/ca_right"]), read(fid["equil/psio"]), read(fid["info/mn_index"]) end @@ -311,11 +319,8 @@ function plot_delta_prime(h5path; save_path=nothing) chi1 = 2π * psio dp_vals = ComplexF64[] - labels = String[] - for s in 1:msing q_s = q_sing[s] - # Find the resonant mode index: (m, n) in mn_index with m/n ≈ q_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 @@ -326,35 +331,36 @@ function plot_delta_prime(h5path; save_path=nothing) dp = (ca_r[resnum, resnum, 2, s] - ca_l[resnum, resnum, 2, s]) / (2π * chi1) push!(dp_vals, dp) end - push!(labels, "q=$(round(q_s, digits=3))") end - colors = [real(v) > 0 ? :red : :steelblue for v in dp_vals] + dp_real = real.(dp_vals) + colors = [v > 0 ? :red : :steelblue for v in dp_real] - p = bar( - 1:msing, abs.(dp_vals); - xticks=(1:msing, labels), - xlabel="rational surface", - ylabel="|Δ'|", - title="Tearing stability Δ' (FFS asymptotic coefficients)", + p = scatter( + psi_sing, dp_real; + xlabel=L"\psi_N", + ylabel="Re(Δ')", + title="Tearing stability Δ'", legend=false, color=colors, - xrotation=30, - bottom_margin=10Plots.mm, - left_margin=5Plots.mm + markersize=7, + markerstrokewidth=0, + left_margin=5Plots.mm, + bottom_margin=5Plots.mm ) + hline!(p, [0.0]; linestyle=:dash, color=:black, label=nothing) isnothing(save_path) || savefig(p, save_path) return p end """ - plot_dcon_summary(h5path; save_path=nothing) + plot_ffs_summary(h5path; save_path=nothing) Four-panel summary of ForceFreeStates (DCON-style) stability results, combining: - Energy eigenvector heatmap (`plot_energy_eigenvectors`) - - Stability criterion vs ψ_N (`plot_stability_criterion`) + - Fixed-boundary stability criterion |D_c| vs ψ_N (`plot_stability_criterion`) - Eigenvalue spectrum (`plot_eigenvalue_spectrum`) - Tearing stability Δ' at each rational surface (`plot_delta_prime`) @@ -373,7 +379,7 @@ panels are shown. A `Plots.jl` plot object. """ -function plot_dcon_summary(h5path; save_path=nothing) +function plot_ffs_summary(h5path; save_path=nothing) has_vac = h5open(h5path, "r") do fid haskey(fid, "vacuum/wt") && !isempty(read(fid["vacuum/wt"])) end @@ -383,8 +389,8 @@ function plot_dcon_summary(h5path; save_path=nothing) if has_vac p_evec = plot_energy_eigenvectors(h5path; matrix_type=:total) - p_eval = plot_eigenvalue_spectrum(h5path; matrix_type=:total) - p = plot(p_evec, p_crit, p_eval, p_dp; layout=(2, 2), size=(1100, 900)) + p_modes = plot_mode_displacement(h5path) + p = plot(p_evec, p_crit, p_modes, p_dp; layout=(2, 2), size=(1100, 900)) 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)) @@ -400,7 +406,7 @@ end Two-panel summary of singular (rational) surface locations and tearing stability: - Left: q(ψ) profile with vertical markers at each rational surface - - Right: `|Δ'|` bar chart at each rational surface (`plot_delta_prime`) + - Right: `Re(Δ')` scatter plot at each rational surface vs ψ_N (`plot_delta_prime`) ### Arguments @@ -423,7 +429,7 @@ function plot_singular_surfaces(h5path; save_path=nothing) p_q = plot( xs, q_profile; - xlabel="ψ_N", + xlabel=L"\psi_N", ylabel="q", title="Safety factor and rational surfaces", legend=false, diff --git a/src/Analysis/PerturbedEquilibrium.jl b/src/Analysis/PerturbedEquilibrium.jl index 5f5b24468..158b9509e 100644 --- a/src/Analysis/PerturbedEquilibrium.jl +++ b/src/Analysis/PerturbedEquilibrium.jl @@ -7,6 +7,7 @@ in the `perturbed_equilibrium/` group of a GPEC HDF5 output file. module PerturbedEquilibrium using HDF5 +using LaTeXStrings using Plots # Check that a PE dataset exists and is non-empty. @@ -19,8 +20,8 @@ end """ plot_resonant_flux(h5path; save_path=nothing) -Bar chart of `|Φ_res|` (normalized resonant flux) per singular surface. One bar series per -toroidal mode n, labeled by n value. +Scatter plot of `|Φ_res|` (normalized resonant flux) per singular surface vs ψ_N. One +marker series per toroidal mode n. Integer-valued q rational surfaces are annotated. Requires the perturbed equilibrium module to have been run and `singular_coupling/resonant_flux` to be present in the HDF5 file. @@ -42,27 +43,28 @@ function plot_resonant_flux(h5path; save_path=nothing) _has_pe_data(h5path, key) || return plot(; title="No resonant flux data — run with perturbed equilibrium enabled", legend=false) - resonant_flux, q_sing, msing, pe_n = h5open(h5path, "r") do fid - read(fid[key]), read(fid["singular/q"]), + resonant_flux, psi_sing, q_sing, msing, pe_n = h5open(h5path, "r") do fid + read(fid[key]), read(fid["singular/psi"]), read(fid["singular/q"]), read(fid["singular/msing"]), read(fid["perturbed_equilibrium/forcing_modes/n"]) end - labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] - xticks_arg = (1:msing, labels) - - p = plot(; xlabel="rational surface", ylabel="|Φ_res|", + p = plot(; xlabel=L"\psi_N", ylabel="|Φ_res|", title="Resonant flux |Φ_res| per surface", legend=:outertopright, - left_margin=5Plots.mm, bottom_margin=10Plots.mm) + left_margin=5Plots.mm, bottom_margin=5Plots.mm) n_vals = unique(pe_n) for nn in n_vals n_rows = findall(==(nn), pe_n) - # Sum over poloidal modes for this n rf_n = [sum(abs.(resonant_flux[n_rows, s])) for s in 1:msing] - bar!(p, 1:msing, rf_n; label="n=$nn", alpha=0.7) + scatter!(p, psi_sing, rf_n; label="n=$nn", markersize=7, markerstrokewidth=0) + # Annotate integer-q surfaces only (avoids overcrowding for n>1) + for s in 1:msing + abs(q_sing[s] - round(q_sing[s])) < 0.05 || continue + annotate!(p, psi_sing[s], rf_n[s], + text(" q=$(round(Int, q_sing[s]))", 8, :left, :black)) + end end - plot!(p; xticks=xticks_arg, xrotation=30) isnothing(save_path) || savefig(p, save_path) return p @@ -71,7 +73,8 @@ end """ plot_island_widths(h5path; save_path=nothing) -Bar chart of island half-width `w/2` per singular surface. +Scatter plot of island half-width `w/2` per singular surface vs ψ_N. +Integer-valued q rational surfaces are annotated. Requires `singular_coupling/island_half_width` in the HDF5 file. @@ -92,24 +95,28 @@ function plot_island_widths(h5path; save_path=nothing) _has_pe_data(h5path, key) || return plot(; title="No island width data — run with perturbed equilibrium enabled", legend=false) - island_hw, q_sing, msing = h5open(h5path, "r") do fid - read(fid[key]), read(fid["singular/q"]), read(fid["singular/msing"]) + island_hw, psi_sing, q_sing, msing = h5open(h5path, "r") do fid + read(fid[key]), read(fid["singular/psi"]), read(fid["singular/q"]), + read(fid["singular/msing"]) end - labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] - - p = bar( - 1:msing, island_hw; - xticks=(1:msing, labels), - xlabel="rational surface", + p = scatter( + psi_sing, island_hw; + xlabel=L"\psi_N", ylabel="w/2", title="Island half-widths", legend=false, color=:steelblue, - xrotation=30, + markersize=7, + markerstrokewidth=0, left_margin=5Plots.mm, - bottom_margin=10Plots.mm + bottom_margin=5Plots.mm ) + for s in 1:msing + abs(q_sing[s] - round(q_sing[s])) < 0.05 || continue + annotate!(p, psi_sing[s], island_hw[s], + text(" q=$(round(Int, q_sing[s]))", 8, :left, :black)) + end isnothing(save_path) || savefig(p, save_path) return p @@ -118,8 +125,9 @@ end """ plot_chirikov_parameter(h5path; save_path=nothing) -Bar chart of the Chirikov overlap parameter per singular surface, with a horizontal reference -line at K = 1 (island overlap threshold). Bars are colored red when K > 1. +Scatter plot of the Chirikov overlap parameter per singular 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 `singular_coupling/chirikov_parameter` in the HDF5 file. @@ -140,26 +148,31 @@ function plot_chirikov_parameter(h5path; save_path=nothing) _has_pe_data(h5path, key) || return plot(; title="No Chirikov data — run with perturbed equilibrium enabled", legend=false) - chirikov, q_sing, msing = h5open(h5path, "r") do fid - read(fid[key]), read(fid["singular/q"]), read(fid["singular/msing"]) + chirikov, psi_sing, q_sing, msing = h5open(h5path, "r") do fid + read(fid[key]), read(fid["singular/psi"]), read(fid["singular/q"]), + read(fid["singular/msing"]) end - labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] colors = [k > 1.0 ? :red : :steelblue for k in chirikov] - p = bar( - 1:msing, chirikov; - xticks=(1:msing, labels), - xlabel="rational surface", + p = scatter( + psi_sing, chirikov; + xlabel=L"\psi_N", ylabel="K (Chirikov)", title="Chirikov overlap parameter", legend=false, color=colors, - xrotation=30, + markersize=7, + markerstrokewidth=0, left_margin=5Plots.mm, - bottom_margin=10Plots.mm + bottom_margin=5Plots.mm ) hline!(p, [1.0]; linestyle=:dash, color=:black, label=nothing) + for s in 1:msing + abs(q_sing[s] - round(q_sing[s])) < 0.05 || continue + annotate!(p, psi_sing[s], chirikov[s], + text(" q=$(round(Int, q_sing[s]))", 8, :left, :black)) + end isnothing(save_path) || savefig(p, save_path) return p @@ -168,8 +181,9 @@ end """ plot_pe_delta_prime(h5path; save_path=nothing) -Bar chart of `|Δ'|` per singular surface computed by the perturbed equilibrium module -(from `singular_coupling/delta_prime`). One bar series per toroidal mode n. +Scatter plot of `Re(Δ')` per singular surface vs ψ_N, computed by the perturbed equilibrium +module (from `singular_coupling/delta_prime`). One marker series per toroidal mode n. +Integer-valued q rational surfaces are annotated. This is complementary to `Analysis.ForceFreeStates.plot_delta_prime`, which uses the FFS asymptotic coefficients. The PE result includes the vacuum Green's function contribution. @@ -193,24 +207,30 @@ function plot_pe_delta_prime(h5path; save_path=nothing) _has_pe_data(h5path, key) || return plot(; title="No PE Δ' data — run with perturbed equilibrium enabled", legend=false) - delta_prime, q_sing, msing, pe_n = h5open(h5path, "r") do fid - read(fid[key]), read(fid["singular/q"]), read(fid["singular/msing"]), + delta_prime, psi_sing, q_sing, msing, pe_n = h5open(h5path, "r") do fid + read(fid[key]), read(fid["singular/psi"]), read(fid["singular/q"]), + read(fid["singular/msing"]), read(fid["perturbed_equilibrium/forcing_modes/n"]) end - labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] - - p = plot(; xlabel="rational surface", ylabel="|Δ'|", + p = plot(; xlabel=L"\psi_N", ylabel="Re(Δ')", title="Tearing stability Δ' (PE)", legend=:outertopright, - left_margin=5Plots.mm, bottom_margin=10Plots.mm) + left_margin=5Plots.mm, bottom_margin=5Plots.mm) + hline!(p, [0.0]; linestyle=:dash, color=:black, label=nothing) n_vals = unique(pe_n) for nn in n_vals n_rows = findall(==(nn), pe_n) - dp_n = [maximum(abs.(delta_prime[n_rows, s])) for s in 1:msing] - bar!(p, 1:msing, dp_n; label="n=$nn", alpha=0.7) + dp_n = [real(delta_prime[n_rows[1], s]) for s in 1:msing] + colors = [v > 0 ? :red : :steelblue for v in dp_n] + scatter!(p, psi_sing, dp_n; label="n=$nn", color=colors, + markersize=7, markerstrokewidth=0) + for s in 1:msing + abs(q_sing[s] - round(q_sing[s])) < 0.05 || continue + annotate!(p, psi_sing[s], dp_n[s], + text(" q=$(round(Int, q_sing[s]))", 8, :left, :black)) + end end - plot!(p; xticks=(1:msing, labels), xrotation=30) isnothing(save_path) || savefig(p, save_path) return p @@ -219,10 +239,10 @@ end """ plot_resonant_field(h5path; save_path=nothing) -Five-panel summary of resonant coupling quantities at each singular surface: +Five-panel summary of resonant coupling quantities at each singular surface vs ψ_N: - `|Φ_res|`: resonant flux (`plot_resonant_flux`) - - `|Δ'|`: tearing stability parameter (`plot_pe_delta_prime`) + - `Re(Δ')`: tearing stability parameter (`plot_pe_delta_prime`) - `|I_res|`: resonant current - `w/2`: island half-width (`plot_island_widths`) - `K`: Chirikov overlap parameter (`plot_chirikov_parameter`) @@ -254,30 +274,33 @@ function plot_resonant_field(h5path; save_path=nothing) return p end -# Internal helper — resonant current bar chart +# Internal helper — resonant current scatter plot function _plot_resonant_current(h5path) key = "perturbed_equilibrium/singular_coupling/resonant_current" _has_pe_data(h5path, key) || return plot(; title="No resonant current data", legend=false) - resonant_current, q_sing, msing, pe_n = h5open(h5path, "r") do fid - read(fid[key]), read(fid["singular/q"]), read(fid["singular/msing"]), + resonant_current, psi_sing, q_sing, msing, pe_n = h5open(h5path, "r") do fid + read(fid[key]), read(fid["singular/psi"]), read(fid["singular/q"]), + read(fid["singular/msing"]), read(fid["perturbed_equilibrium/forcing_modes/n"]) end - labels = ["q=$(round(q_sing[s], digits=3))" for s in 1:msing] - - p = plot(; xlabel="rational surface", ylabel="|I_res|", + p = plot(; xlabel=L"\psi_N", ylabel="|I_res|", title="Resonant current |I_res| per surface", legend=:outertopright, - left_margin=5Plots.mm, bottom_margin=10Plots.mm) + left_margin=5Plots.mm, bottom_margin=5Plots.mm) n_vals = unique(pe_n) for nn in n_vals n_rows = findall(==(nn), pe_n) rc_n = [sum(abs.(resonant_current[n_rows, s])) for s in 1:msing] - bar!(p, 1:msing, rc_n; label="n=$nn", alpha=0.7) + scatter!(p, psi_sing, rc_n; label="n=$nn", markersize=7, markerstrokewidth=0) + for s in 1:msing + abs(q_sing[s] - round(q_sing[s])) < 0.05 || continue + annotate!(p, psi_sing[s], rc_n[s], + text(" q=$(round(Int, q_sing[s]))", 8, :left, :black)) + end end - plot!(p; xticks=(1:msing, labels), xrotation=30) return p end @@ -342,7 +365,7 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) # Top panel: line plot per mode — only label resonant range m ∈ [0, nhigh·q95) m_max_legend = nhigh * q95 p1 = plot(; - xlabel="ψ_N", + xlabel=L"\psi_N", ylabel="|$(component)|", title="Mode spectrogram: $(component)", legend=:outertopright, @@ -361,7 +384,7 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) p2 = heatmap( collect(m_vals), psi_response, abs.(data_mn); xlabel="m", - ylabel="ψ_N", + ylabel=L"\psi_N", title="", colorbar_title="|$(component)|", left_margin=5Plots.mm, From 5bafd6e4c07b6f9551a7b91e172a7e80df504004 Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Wed, 25 Mar 2026 19:43:40 -0400 Subject: [PATCH 06/12] ANALYSIS - IMPROVEMENT - Fix pe_resonant_field shared x-axis and legend alignment Remove titles from all 5 subplots in plot_resonant_field (ylabels carry the same information). Move outertopright legends inside (:topright) for the top 3 panels so all 5 panels share the same right margin and x-axes align. Suppress x-axis labels and ticks on the top 4 panels; only the bottom (Chirikov) panel shows the x-axis. Co-Authored-By: Claude Sonnet 4.6 --- src/Analysis/Equilibrium.jl | 16 +++++++++------- src/Analysis/ForceFreeStates.jl | 14 +++++++------- src/Analysis/PerturbedEquilibrium.jl | 27 ++++++++++++++++++++------- 3 files changed, 36 insertions(+), 21 deletions(-) diff --git a/src/Analysis/Equilibrium.jl b/src/Analysis/Equilibrium.jl index 46cf3b448..93dd07827 100644 --- a/src/Analysis/Equilibrium.jl +++ b/src/Analysis/Equilibrium.jl @@ -95,7 +95,7 @@ function plot_qprofile(h5path; show_singular=true, save_path=nothing) p = plot( xs, q; - xlabel=L"\psi_N", + xlabel="ψₙ", ylabel="q", title="", legend=false, @@ -148,7 +148,7 @@ function plot_pressure_profile(h5path; save_path=nothing) p = plot( xs, mu0p; - xlabel=L"\psi_N", + xlabel="ψₙ", ylabel="μ₀p", title="", legend=false, @@ -189,7 +189,7 @@ function plot_f_profile(h5path; save_path=nothing) p = plot( xs, twopif; - xlabel=L"\psi_N", + xlabel="ψₙ", ylabel="2πF", title="", legend=false, @@ -274,7 +274,7 @@ function plot_flux_surfaces_h5(h5path; n_psi=11, n_theta=18, save_path=nothing) plot!(p, R_grid[:, itheta], Z_grid[:, itheta]; color=:tomato, linewidth=0.8, label=label) end - # Draw rational surface flux contours in red, annotated at the outboard midplane + # Draw rational surface flux contours in red, annotated at the bottom (min Z) of each for s in 1:msing idx = argmin(abs.(xs_rz .- psi_sing[s])) q_label = abs(q_sing[s] - round(q_sing[s])) < 0.05 ? @@ -282,7 +282,9 @@ function plot_flux_surfaces_h5(h5path; n_psi=11, n_theta=18, save_path=nothing) plot!(p, [R_grid[idx, :]; R_grid[idx, 1]], [Z_grid[idx, :]; Z_grid[idx, 1]]; color=:red, linewidth=1.5, label=s == 1 ? "Rational surface" : "") - annotate!(p, R_grid[idx, 1], Z_grid[idx, 1], text(" $q_label", 7, :left, :red)) + itheta_bot = argmin(Z_grid[idx, :]) + 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) @@ -322,7 +324,7 @@ function plot_gse_by_theta(h5path; n_theta_lines=8, save_path=nothing) theta_indices = round.(Int, range(1, ntheta; length=n_theta_lines)) p = plot(; - xlabel=L"\psi_N", + xlabel="ψₙ", ylabel="GSE error", title="", yscale=:log10, @@ -379,7 +381,7 @@ function plot_gse_integrated(h5path; save_path=nothing) p = plot( xs, vec(errlogi); - xlabel=L"\psi_N", + xlabel="ψₙ", ylabel="log₁₀(integrated GSE)", title="Flux-surface-integrated Grad-Shafranov error", legend=false, diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index 08683595a..c8422bd7a 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -39,11 +39,11 @@ function plot_mode_displacement(h5path; modes=1:5, save_path=nothing) mhigh = mlow + mpert - 1 dW = isempty(et) ? nothing : et[1] title_str = isnothing(dW) ? "Least stable mode" : - "Least stable mode, dW = $(round(real(dW), sigdigits=4))" + "Least stable mode, δW = $(round(real(dW), sigdigits=4))" p = plot(; - xlabel=L"\psi_N", - ylabel=L"|\xi_\psi|", + xlabel="ψₙ", + ylabel="|ξ^ψ|", title=title_str, left_margin=5Plots.mm, bottom_margin=5Plots.mm @@ -166,8 +166,8 @@ function plot_stability_criterion(h5path; save_path=nothing) p = plot( psi, crit; - xlabel=L"\psi_N", - ylabel=L"|D_c|", + xlabel="ψₙ", + ylabel="|Dᶜ|", title="Fixed-boundary stability", legend=false, left_margin=5Plots.mm, @@ -338,7 +338,7 @@ function plot_delta_prime(h5path; save_path=nothing) p = scatter( psi_sing, dp_real; - xlabel=L"\psi_N", + xlabel="ψₙ", ylabel="Re(Δ')", title="Tearing stability Δ'", legend=false, @@ -429,7 +429,7 @@ function plot_singular_surfaces(h5path; save_path=nothing) p_q = plot( xs, q_profile; - xlabel=L"\psi_N", + xlabel="ψₙ", ylabel="q", title="Safety factor and rational surfaces", legend=false, diff --git a/src/Analysis/PerturbedEquilibrium.jl b/src/Analysis/PerturbedEquilibrium.jl index 158b9509e..298e2d1c1 100644 --- a/src/Analysis/PerturbedEquilibrium.jl +++ b/src/Analysis/PerturbedEquilibrium.jl @@ -49,7 +49,7 @@ function plot_resonant_flux(h5path; save_path=nothing) read(fid["perturbed_equilibrium/forcing_modes/n"]) end - p = plot(; xlabel=L"\psi_N", ylabel="|Φ_res|", + p = plot(; xlabel="ψₙ", ylabel="|Φ_res|", title="Resonant flux |Φ_res| per surface", legend=:outertopright, left_margin=5Plots.mm, bottom_margin=5Plots.mm) @@ -102,7 +102,7 @@ function plot_island_widths(h5path; save_path=nothing) p = scatter( psi_sing, island_hw; - xlabel=L"\psi_N", + xlabel="ψₙ", ylabel="w/2", title="Island half-widths", legend=false, @@ -157,7 +157,7 @@ function plot_chirikov_parameter(h5path; save_path=nothing) p = scatter( psi_sing, chirikov; - xlabel=L"\psi_N", + xlabel="ψₙ", ylabel="K (Chirikov)", title="Chirikov overlap parameter", legend=false, @@ -213,7 +213,7 @@ function plot_pe_delta_prime(h5path; save_path=nothing) read(fid["perturbed_equilibrium/forcing_modes/n"]) end - p = plot(; xlabel=L"\psi_N", ylabel="Re(Δ')", + p = plot(; xlabel="ψₙ", ylabel="Re(Δ')", title="Tearing stability Δ' (PE)", legend=:outertopright, left_margin=5Plots.mm, bottom_margin=5Plots.mm) hline!(p, [0.0]; linestyle=:dash, color=:black, label=nothing) @@ -268,6 +268,19 @@ function plot_resonant_field(h5path; save_path=nothing) p4 = plot_island_widths(h5path) p5 = plot_chirikov_parameter(h5path) + # Move outer legends inside so all panels have the same right margin + for p in (p1, p2, p3) + plot!(p; title="", legend=:topright) + end + for p in (p4, p5) + plot!(p; title="") + end + + hide_xaxis!(p) = plot!(p; xlabel="", xformatter=_->"", bottom_margin=1Plots.mm) + for p in (p1, p2, p3, p4) + hide_xaxis!(p) + end + p = plot(p1, p2, p3, p4, p5; layout=(5, 1), size=(900, 1600)) isnothing(save_path) || savefig(p, save_path) @@ -286,7 +299,7 @@ function _plot_resonant_current(h5path) read(fid["perturbed_equilibrium/forcing_modes/n"]) end - p = plot(; xlabel=L"\psi_N", ylabel="|I_res|", + p = plot(; xlabel="ψₙ", ylabel="|I_res|", title="Resonant current |I_res| per surface", legend=:outertopright, left_margin=5Plots.mm, bottom_margin=5Plots.mm) @@ -365,7 +378,7 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) # Top panel: line plot per mode — only label resonant range m ∈ [0, nhigh·q95) m_max_legend = nhigh * q95 p1 = plot(; - xlabel=L"\psi_N", + xlabel="ψₙ", ylabel="|$(component)|", title="Mode spectrogram: $(component)", legend=:outertopright, @@ -384,7 +397,7 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) p2 = heatmap( collect(m_vals), psi_response, abs.(data_mn); xlabel="m", - ylabel=L"\psi_N", + ylabel="ψₙ", title="", colorbar_title="|$(component)|", left_margin=5Plots.mm, From 065d7f1a4cc31979e768a4bfbdb958a5aa9079eb Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Thu, 26 Mar 2026 16:38:44 -0400 Subject: [PATCH 07/12] Analyze - WIP - Cleans up a bunch of extras claude added for no good reason and fixes some labels, xlims, etc. --- .../analyze_example.jl | 14 +- .../DIIID-like_ideal_example/make_plots.jl | 73 -------- .../Solovev_ideal_example/analyze_example.jl | 9 +- .../analyze_example.jl | 14 +- src/Analysis/Equilibrium.jl | 75 +------- src/Analysis/ForceFreeStates.jl | 163 ++---------------- src/Analysis/PerturbedEquilibrium.jl | 22 +-- 7 files changed, 46 insertions(+), 324 deletions(-) delete mode 100644 examples/DIIID-like_ideal_example/make_plots.jl diff --git a/examples/DIIID-like_ideal_example/analyze_example.jl b/examples/DIIID-like_ideal_example/analyze_example.jl index 7a562499b..dbf5da961 100644 --- a/examples/DIIID-like_ideal_example/analyze_example.jl +++ b/examples/DIIID-like_ideal_example/analyze_example.jl @@ -5,16 +5,8 @@ using GeneralizedPerturbedEquilibrium: Analysis plotlyjs() h5path = joinpath(@__DIR__, "gpec.h5") -p_modes = Analysis.ForceFreeStates.plot_mode_displacement(h5path; modes=1:5) -p_eigen = Analysis.ForceFreeStates.plot_eigenmode_summary(h5path) -p_stab = Analysis.ForceFreeStates.plot_stability_criterion(h5path) # Summary plots -p_dcon_summary = Analysis.ForceFreeStates.plot_dcon_summary(h5path) -p_sing = Analysis.ForceFreeStates.plot_singular_surfaces(h5path) -p_equil = Analysis.Equilibrium.plot_equilibrium_summary(h5path) - -# Perturbed equilibrium plots (requires PE module run with forcing.dat) -p_resfield = Analysis.PerturbedEquilibrium.plot_resonant_field(h5path) -p_pe_summary = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path) -p_spectrogram = Analysis.PerturbedEquilibrium.plot_mode_spectrogram(h5path; component=:xi_psi) +p_eq = Analysis.Equilibrium.plot_equilibrium_summary(h5path) +p_ffs = Analysis.ForceFreeStates.plot_ffs_summary(h5path) +p_pe = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path) \ No newline at end of file diff --git a/examples/DIIID-like_ideal_example/make_plots.jl b/examples/DIIID-like_ideal_example/make_plots.jl deleted file mode 100644 index 68e728e0a..000000000 --- a/examples/DIIID-like_ideal_example/make_plots.jl +++ /dev/null @@ -1,73 +0,0 @@ -using Pkg -Pkg.activate(joinpath(@__DIR__, "../..")) -using GeneralizedPerturbedEquilibrium, Plots, GeneralizedPerturbedEquilibrium.Analysis -gr() - -h5path = joinpath(@__DIR__, "gpec.h5") -outdir = joinpath(@__DIR__, "plots") -mkpath(outdir) - -# --- Equilibrium --- -p = Analysis.Equilibrium.plot_equilibrium_summary(h5path) -savefig(p, joinpath(outdir, "equil_summary.png")) - -p = Analysis.Equilibrium.plot_flux_surfaces_h5(h5path) -savefig(p, joinpath(outdir, "equil_flux_surfaces.png")) - -p = Analysis.Equilibrium.plot_qprofile(h5path) -savefig(p, joinpath(outdir, "equil_qprofile.png")) - -p = Analysis.Equilibrium.plot_pressure_profile(h5path) -savefig(p, joinpath(outdir, "equil_pressure.png")) - -p = Analysis.Equilibrium.plot_f_profile(h5path) -savefig(p, joinpath(outdir, "equil_F.png")) - -# --- ForceFreeStates --- -p = Analysis.ForceFreeStates.plot_mode_displacement(h5path; modes=1:5) -savefig(p, joinpath(outdir, "ffs_mode_displacement.png")) - -p = Analysis.ForceFreeStates.plot_stability_criterion(h5path) -savefig(p, joinpath(outdir, "ffs_stability_criterion.png")) - -p = Analysis.ForceFreeStates.plot_energy_eigenvectors(h5path) -savefig(p, joinpath(outdir, "ffs_energy_eigenvectors.png")) - -p = Analysis.ForceFreeStates.plot_eigenvalue_spectrum(h5path) -savefig(p, joinpath(outdir, "ffs_eigenvalue_spectrum.png")) - -p = Analysis.ForceFreeStates.plot_delta_prime(h5path) -savefig(p, joinpath(outdir, "ffs_delta_prime.png")) - -p = Analysis.ForceFreeStates.plot_ffs_summary(h5path) -savefig(p, joinpath(outdir, "ffs_summary.png")) - -p = Analysis.ForceFreeStates.plot_singular_surfaces(h5path) -savefig(p, joinpath(outdir, "ffs_singular_surfaces.png")) - -# --- PerturbedEquilibrium --- -p = Analysis.PerturbedEquilibrium.plot_resonant_flux(h5path) -savefig(p, joinpath(outdir, "pe_resonant_flux.png")) - -p = Analysis.PerturbedEquilibrium.plot_island_widths(h5path) -savefig(p, joinpath(outdir, "pe_island_widths.png")) - -p = Analysis.PerturbedEquilibrium.plot_chirikov_parameter(h5path) -savefig(p, joinpath(outdir, "pe_chirikov.png")) - -p = Analysis.PerturbedEquilibrium.plot_pe_delta_prime(h5path) -savefig(p, joinpath(outdir, "pe_delta_prime.png")) - -p = Analysis.PerturbedEquilibrium.plot_resonant_field(h5path) -savefig(p, joinpath(outdir, "pe_resonant_field.png")) - -p = Analysis.PerturbedEquilibrium.plot_mode_spectrogram(h5path; component=:xi_psi) -savefig(p, joinpath(outdir, "pe_spectrogram_xi_psi.png")) - -p = Analysis.PerturbedEquilibrium.plot_mode_spectrogram(h5path; component=:b_psi) -savefig(p, joinpath(outdir, "pe_spectrogram_b_psi.png")) - -p = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path) -savefig(p, joinpath(outdir, "pe_summary.png")) - -println("All plots saved to: $outdir") diff --git a/examples/Solovev_ideal_example/analyze_example.jl b/examples/Solovev_ideal_example/analyze_example.jl index 45a2754e4..dbf5da961 100644 --- a/examples/Solovev_ideal_example/analyze_example.jl +++ b/examples/Solovev_ideal_example/analyze_example.jl @@ -5,11 +5,8 @@ using GeneralizedPerturbedEquilibrium: Analysis plotlyjs() h5path = joinpath(@__DIR__, "gpec.h5") -p_modes = Analysis.ForceFreeStates.plot_mode_displacement(h5path; modes=1:5) -p_eigen = Analysis.ForceFreeStates.plot_eigenmode_summary(h5path) -p_stab = Analysis.ForceFreeStates.plot_stability_criterion(h5path) # Summary plots -p_dcon_summary = Analysis.ForceFreeStates.plot_dcon_summary(h5path) -p_sing = Analysis.ForceFreeStates.plot_singular_surfaces(h5path) -p_equil = Analysis.Equilibrium.plot_equilibrium_summary(h5path) +p_eq = Analysis.Equilibrium.plot_equilibrium_summary(h5path) +p_ffs = Analysis.ForceFreeStates.plot_ffs_summary(h5path) +p_pe = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path) \ No newline at end of file diff --git a/examples/Solovev_ideal_example_multi_n/analyze_example.jl b/examples/Solovev_ideal_example_multi_n/analyze_example.jl index b4e4e2539..ec1873cec 100644 --- a/examples/Solovev_ideal_example_multi_n/analyze_example.jl +++ b/examples/Solovev_ideal_example_multi_n/analyze_example.jl @@ -6,17 +6,13 @@ plotlyjs() # Analyze the multi-n run h5path = joinpath(@__DIR__, "gpec.h5") -p_modes = Analysis.ForceFreeStates.plot_mode_displacement(h5path; modes=1:5) -p_eigen = Analysis.ForceFreeStates.plot_eigenmode_summary(h5path) -p_stab = Analysis.ForceFreeStates.plot_stability_criterion(h5path) - # Summary plots -p_dcon_summary = Analysis.ForceFreeStates.plot_dcon_summary(h5path) -p_sing = Analysis.ForceFreeStates.plot_singular_surfaces(h5path) -p_equil = Analysis.Equilibrium.plot_equilibrium_summary(h5path) +p_eq = Analysis.Equilibrium.plot_equilibrium_summary(h5path) +p_ffs = Analysis.ForceFreeStates.plot_ffs_summary(h5path) +p_pe = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path) # Analyze the single-n runs h5path_n1 = joinpath(@__DIR__, "single_n_1", "euler_n1.h5") h5path_n2 = joinpath(@__DIR__, "single_n_2", "euler_n2.h5") -p_modes_n1 = Analysis.ForceFreeStates.plot_mode_displacement(h5path_n1; modes=1:5) -p_modes_n2 = Analysis.ForceFreeStates.plot_mode_displacement(h5path_n2; modes=1:5) +p_n1 = Analysis.ForceFreeStates.plot_ffs_summary(h5path_n1) +p_n2 = Analysis.ForceFreeStates.plot_ffs_summary(h5path_n2) diff --git a/src/Analysis/Equilibrium.jl b/src/Analysis/Equilibrium.jl index 93dd07827..66a2466e8 100644 --- a/src/Analysis/Equilibrium.jl +++ b/src/Analysis/Equilibrium.jl @@ -9,65 +9,6 @@ using HDF5 using LaTeXStrings using Plots -""" - plot_flux_surfaces(plasma_eq; n_psi=11, n_theta=13) - -Plot flux surface contours (constant ψ, blue) and field-line angle contours (constant θ, red) -in physical (R, Z) space. - -### Arguments - - - `plasma_eq`: A `GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium` object - -### Keyword arguments - - - `n_psi`: Number of constant-ψ contours to draw (default: 11) - - `n_theta`: Number of constant-θ contours to draw (default: 13) - -### Returns - -A `Plots.jl` plot object. -""" -function plot_flux_surfaces(plasma_eq; n_psi=11, n_theta=13) - n_psi_grid = size(plasma_eq.rzphi_rsquared.nodal_derivs.partials, 2) - n_theta_grid = size(plasma_eq.rzphi_rsquared.nodal_derivs.partials, 3) - - # Build R and 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 - rfac = @. sqrt(max(0.0, plasma_eq.rzphi_rsquared.nodal_derivs.partials[1, ipsi, :])) - angle = @. 2π * (plasma_eq.rzphi_ys + plasma_eq.rzphi_offset.nodal_derivs.partials[1, ipsi, :]) - R_grid[ipsi, :] = plasma_eq.ro .+ rfac .* cos.(angle) - Z_grid[ipsi, :] = plasma_eq.zo .+ rfac .* sin.(angle) - end - - p = plot(; - title="Flux Coordinate System Contours in (R, Z)", - xlabel="R [m]", - ylabel="Z [m]", - aspect_ratio=:equal, - legend=:outertopright - ) - - psi_indices = round.(Int, range(1, n_psi_grid; length=n_psi)) - theta_indices = round.(Int, range(1, n_theta_grid; length=n_theta)) - - for (i, ipsi) in enumerate(psi_indices) - label = i == 1 ? "Constant ψ" : "" - plot!(p, [R_grid[ipsi, :]; R_grid[ipsi, 1]], [Z_grid[ipsi, :]; Z_grid[ipsi, 1]]; - color=:blue, linewidth=1.5, label=label) - end - - for (i, itheta) in enumerate(theta_indices) - label = i == 1 ? "Constant θ" : "" - plot!(p, R_grid[:, itheta], Z_grid[:, itheta]; - color=:red, linewidth=1.0, label=label) - end - - return p -end - """ plot_qprofile(h5path; show_singular=true, save_path=nothing) @@ -95,7 +36,7 @@ function plot_qprofile(h5path; show_singular=true, save_path=nothing) p = plot( xs, q; - xlabel="ψₙ", + xlabel="Norm. Poloidal Flux", ylabel="q", title="", legend=false, @@ -148,7 +89,7 @@ function plot_pressure_profile(h5path; save_path=nothing) p = plot( xs, mu0p; - xlabel="ψₙ", + xlabel="Norm. Poloidal Flux", ylabel="μ₀p", title="", legend=false, @@ -189,7 +130,7 @@ function plot_f_profile(h5path; save_path=nothing) p = plot( xs, twopif; - xlabel="ψₙ", + xlabel="Norm. Poloidal Flux", ylabel="2πF", title="", legend=false, @@ -206,7 +147,7 @@ function plot_f_profile(h5path; save_path=nothing) end """ - plot_flux_surfaces_h5(h5path; n_psi=11, n_theta=18, save_path=nothing) + plot_flux_surfaces(h5path; 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. @@ -228,7 +169,7 @@ Theta spokes are drawn at `n_theta` evenly spaced values. A `Plots.jl` plot object. """ -function plot_flux_surfaces_h5(h5path; n_psi=11, n_theta=18, save_path=nothing) +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["splines/rzphi/rcoords"]), read(fid["splines/rzphi/offset"]), read(fid["splines/rzphi/xs"]), read(fid["splines/rzphi/ys"]), @@ -324,7 +265,7 @@ function plot_gse_by_theta(h5path; n_theta_lines=8, save_path=nothing) theta_indices = round.(Int, range(1, ntheta; length=n_theta_lines)) p = plot(; - xlabel="ψₙ", + xlabel="Norm. Poloidal Flux", ylabel="GSE error", title="", yscale=:log10, @@ -381,7 +322,7 @@ function plot_gse_integrated(h5path; save_path=nothing) p = plot( xs, vec(errlogi); - xlabel="ψₙ", + xlabel="Norm. Poloidal Flux", ylabel="log₁₀(integrated GSE)", title="Flux-surface-integrated Grad-Shafranov error", legend=false, @@ -427,7 +368,7 @@ function plot_equilibrium_summary(h5path; save_path=nothing) 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_h5(h5path) + 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) diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index c8422bd7a..610744437 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -42,7 +42,8 @@ function plot_mode_displacement(h5path; modes=1:5, save_path=nothing) "Least stable mode, δW = $(round(real(dW), sigdigits=4))" p = plot(; - xlabel="ψₙ", + xlims=(0, 1), + xlabel="Norm. Poloidal Flux", ylabel="|ξ^ψ|", title=title_str, left_margin=5Plots.mm, @@ -57,92 +58,9 @@ function plot_mode_displacement(h5path; modes=1:5, save_path=nothing) return p end -""" - plot_eigenmode_summary(h5path; save_path=nothing) - -Three-panel summary of the free-boundary energy matrix eigenmodes, analogous to the -DCON summary plot produced by OMFIT GPEC. - -Panels: - - - Top: |eigenvector| vs index for the least stable mode - - Bottom-left: heatmap of |W_t| (eigenvectors) vs mode index - - Bottom-right: |eigenvalue| on log scale vs mode index - -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`) - -### Keyword arguments - - - `save_path`: If provided, save the figure to this path (default: `nothing`) - -### Returns - -A `Plots.jl` plot object. -""" -function plot_eigenmode_summary(h5path; save_path=nothing) - wt, et, psio, mlow = h5open(h5path, "r") do fid - read(fid["vacuum/wt"]), read(fid["vacuum/et"]), - read(fid["equil/psio"]), read(fid["info/mlow"]) - end - - isempty(wt) && error("No vacuum data in $h5path; rerun with vac_flag = true") - - chi1 = 2π * psio - wt = wt * (chi1 * 1e-3) - - nmn = size(wt, 1) - nmodes = size(wt, 2) - m_vals = (0:(nmn-1)) .+ mlow - - p1 = plot( - m_vals, abs.(wt[:, 1]); - xlabel="m", - ylabel="|Eigenvector|", - title="Mode 1, |λ₁| = $(round(abs(et[1]), digits=3))", - legend=false, - left_margin=5Plots.mm, - bottom_margin=5Plots.mm - ) - - p2 = heatmap( - m_vals, 1:nmodes, abs.(wt'); - xlabel="m", - ylabel="mode index", - colorbar_title="|Wₜ|", - left_margin=5Plots.mm, - right_margin=10Plots.mm, - bottom_margin=5Plots.mm - ) - - colors = [real(e) < 0 ? :red : :blue for e in et] - p3 = scatter( - abs.(et), 1:nmodes; - xlabel="|Eigenvalue|", - ylabel="mode index", - xscale=:log10, - legend=false, - color=colors, - markerstrokewidth=0, - left_margin=5Plots.mm, - bottom_margin=5Plots.mm - ) - - # Top panel (p1) shares m-axis with heatmap (p2); blank cell keeps widths aligned. - # Plots.jl assigns plots sequentially to ALL layout slots including `_`, so pass an - # explicit empty plot() to fill the blank top-right cell. - l = @layout [a{0.25h} _{0.25w}; b c{0.25w}] - p = plot(p1, plot(), p2, p3; layout=l, size=(950, 750)) - - isnothing(save_path) || savefig(p, save_path) - return p -end """ - plot_stability_criterion(h5path; save_path=nothing) + plot_fixed_boundarystability_criterion(h5path; save_path=nothing) Plot the stability criterion (smallest eigenvalue of W⁻¹, `crit`) vs ψ_N. A sign change in `crit` during integration indicates an ideal fixed-boundary instability. @@ -159,16 +77,17 @@ A sign change in `crit` during integration indicates an ideal fixed-boundary ins A `Plots.jl` plot object. """ -function plot_stability_criterion(h5path; save_path=nothing) +function plot_fixed_boundarystability_criterion(h5path; save_path=nothing) psi, crit = h5open(h5path, "r") do fid read(fid["integration/psi"]), read(fid["integration/crit"]) end p = plot( psi, crit; - xlabel="ψₙ", + xlims=(0, 1), + xlabel="Norm. Poloidal Flux", ylabel="|Dᶜ|", - title="Fixed-boundary stability", + title="Fixed-Boundary Stability", legend=false, left_margin=5Plots.mm, bottom_margin=5Plots.mm @@ -221,10 +140,10 @@ function plot_energy_eigenvectors(h5path; matrix_type=:total, save_path=nothing) p = heatmap( m_vals, 1:nmodes, abs.(wt_scaled'); - xlabel="m", - ylabel="mode index", - title="Energy eigenvectors |Wₜ| (total)", - colorbar_title="|Wₜ|", + xlabel="Poloidal Harmonic", + ylabel="Eigenmode Index", + title="Total Energy Eigenvectors", + colorbar_title="Harmonic Amplitude", left_margin=5Plots.mm, right_margin=10Plots.mm, bottom_margin=5Plots.mm @@ -235,7 +154,7 @@ function plot_energy_eigenvectors(h5path; matrix_type=:total, save_path=nothing) end """ - plot_eigenvalue_spectrum(h5path; matrix_type=:total, save_path=nothing) + plot_eigenvalues(h5path; 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. @@ -253,7 +172,7 @@ or green (stable, Re < 0), with a dashed reference line at zero. A `Plots.jl` plot object. """ -function plot_eigenvalue_spectrum(h5path; matrix_type=:total, save_path=nothing) +function plot_eigenvalues(h5path; matrix_type=:total, save_path=nothing) dataset = Dict(:total => "vacuum/et", :plasma => "vacuum/ep", :vacuum => "vacuum/ev") haskey(dataset, matrix_type) || error("matrix_type must be :total, :plasma, or :vacuum") @@ -338,7 +257,8 @@ function plot_delta_prime(h5path; save_path=nothing) p = scatter( psi_sing, dp_real; - xlabel="ψₙ", + xlims=(0, 1), + xlabel="Norm. Poloidal Flux", ylabel="Re(Δ')", title="Tearing stability Δ'", legend=false, @@ -361,7 +281,7 @@ Four-panel summary of ForceFreeStates (DCON-style) stability results, combining: - Energy eigenvector heatmap (`plot_energy_eigenvectors`) - Fixed-boundary stability criterion |D_c| vs ψ_N (`plot_stability_criterion`) - - Eigenvalue spectrum (`plot_eigenvalue_spectrum`) + - Eigenvalue spectrum (`plot_eigenvalues`) - Tearing stability Δ' at each rational surface (`plot_delta_prime`) If no vacuum data is present (`vac_flag = false`), only the stability criterion and Δ' @@ -400,55 +320,4 @@ function plot_ffs_summary(h5path; save_path=nothing) return p end -""" - plot_singular_surfaces(h5path; save_path=nothing) - -Two-panel summary of singular (rational) surface locations and tearing stability: - - - Left: q(ψ) profile with vertical markers at each rational surface - - Right: `Re(Δ')` scatter plot at each rational surface vs ψ_N (`plot_delta_prime`) - -### Arguments - - - `h5path`: Path to a GPEC HDF5 output file - -### Keyword arguments - - - `save_path`: If provided, save the figure to this path (default: `nothing`) - -### Returns - -A `Plots.jl` plot object. -""" -function plot_singular_surfaces(h5path; save_path=nothing) - xs, q_profile, q0, q95, msing, psi_sing, q_sing = h5open(h5path, "r") do fid - read(fid["splines/profiles/xs"]), read(fid["splines/profiles/q"]), - read(fid["equil/q0"]), read(fid["equil/q95"]), - read(fid["singular/msing"]), read(fid["singular/psi"]), read(fid["singular/q"]) - end - - p_q = plot( - xs, q_profile; - xlabel="ψₙ", - ylabel="q", - title="Safety factor and rational surfaces", - legend=false, - left_margin=5Plots.mm, - bottom_margin=5Plots.mm - ) - hline!(p_q, [q0, q95]; linestyle=:dot, color=:gray, label=nothing) - for s in 1:msing - vline!(p_q, [psi_sing[s]]; linestyle=:dash, color=:red, label=nothing) - annotate!(p_q, psi_sing[s], q_sing[s], - text(" q=$(round(q_sing[s], digits=2))", 7, :left, :red)) - end - - p_dp = plot_delta_prime(h5path) - - p = plot(p_q, p_dp; layout=(1, 2), size=(1100, 500)) - - isnothing(save_path) || savefig(p, save_path) - return p -end - end # module ForceFreeStates diff --git a/src/Analysis/PerturbedEquilibrium.jl b/src/Analysis/PerturbedEquilibrium.jl index 298e2d1c1..6b2465a26 100644 --- a/src/Analysis/PerturbedEquilibrium.jl +++ b/src/Analysis/PerturbedEquilibrium.jl @@ -49,7 +49,7 @@ function plot_resonant_flux(h5path; save_path=nothing) read(fid["perturbed_equilibrium/forcing_modes/n"]) end - p = plot(; xlabel="ψₙ", ylabel="|Φ_res|", + p = plot(; xlabel="Norm. Poloidal Flux", ylabel="|Φ_res|", title="Resonant flux |Φ_res| per surface", legend=:outertopright, left_margin=5Plots.mm, bottom_margin=5Plots.mm) @@ -102,7 +102,7 @@ function plot_island_widths(h5path; save_path=nothing) p = scatter( psi_sing, island_hw; - xlabel="ψₙ", + xlabel="Norm. Poloidal Flux", ylabel="w/2", title="Island half-widths", legend=false, @@ -157,7 +157,7 @@ function plot_chirikov_parameter(h5path; save_path=nothing) p = scatter( psi_sing, chirikov; - xlabel="ψₙ", + xlabel="Norm. Poloidal Flux", ylabel="K (Chirikov)", title="Chirikov overlap parameter", legend=false, @@ -179,7 +179,7 @@ function plot_chirikov_parameter(h5path; save_path=nothing) end """ - plot_pe_delta_prime(h5path; save_path=nothing) + plot_driven_delta_prime(h5path; save_path=nothing) Scatter plot of `Re(Δ')` per singular surface vs ψ_N, computed by the perturbed equilibrium module (from `singular_coupling/delta_prime`). One marker series per toroidal mode n. @@ -202,7 +202,7 @@ Requires `singular_coupling/delta_prime` in the HDF5 file. A `Plots.jl` plot object. """ -function plot_pe_delta_prime(h5path; save_path=nothing) +function plot_driven_delta_prime(h5path; save_path=nothing) key = "perturbed_equilibrium/singular_coupling/delta_prime" _has_pe_data(h5path, key) || return plot(; title="No PE Δ' data — run with perturbed equilibrium enabled", legend=false) @@ -213,7 +213,7 @@ function plot_pe_delta_prime(h5path; save_path=nothing) read(fid["perturbed_equilibrium/forcing_modes/n"]) end - p = plot(; xlabel="ψₙ", ylabel="Re(Δ')", + p = plot(; xlabel="Norm. Poloidal Flux", ylabel="Re(Δ')", title="Tearing stability Δ' (PE)", legend=:outertopright, left_margin=5Plots.mm, bottom_margin=5Plots.mm) hline!(p, [0.0]; linestyle=:dash, color=:black, label=nothing) @@ -242,7 +242,7 @@ end Five-panel summary of resonant coupling quantities at each singular surface vs ψ_N: - `|Φ_res|`: resonant flux (`plot_resonant_flux`) - - `Re(Δ')`: tearing stability parameter (`plot_pe_delta_prime`) + - `Re(Δ')`: tearing stability parameter (`plot_driven_delta_prime`) - `|I_res|`: resonant current - `w/2`: island half-width (`plot_island_widths`) - `K`: Chirikov overlap parameter (`plot_chirikov_parameter`) @@ -263,7 +263,7 @@ A `Plots.jl` plot object. """ function plot_resonant_field(h5path; save_path=nothing) p1 = plot_resonant_flux(h5path) - p2 = plot_pe_delta_prime(h5path) + p2 = plot_driven_delta_prime(h5path) p3 = _plot_resonant_current(h5path) p4 = plot_island_widths(h5path) p5 = plot_chirikov_parameter(h5path) @@ -299,7 +299,7 @@ function _plot_resonant_current(h5path) read(fid["perturbed_equilibrium/forcing_modes/n"]) end - p = plot(; xlabel="ψₙ", ylabel="|I_res|", + p = plot(; xlabel="Norm. Poloidal Flux", ylabel="|I_res|", title="Resonant current |I_res| per surface", legend=:outertopright, left_margin=5Plots.mm, bottom_margin=5Plots.mm) @@ -378,7 +378,7 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) # Top panel: line plot per mode — only label resonant range m ∈ [0, nhigh·q95) m_max_legend = nhigh * q95 p1 = plot(; - xlabel="ψₙ", + xlabel="Norm. Poloidal Flux", ylabel="|$(component)|", title="Mode spectrogram: $(component)", legend=:outertopright, @@ -397,7 +397,7 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) p2 = heatmap( collect(m_vals), psi_response, abs.(data_mn); xlabel="m", - ylabel="ψₙ", + ylabel="Norm. Poloidal Flux", title="", colorbar_title="|$(component)|", left_margin=5Plots.mm, From 79760c677deac6995f86af139c54ed4c51997a92 Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Thu, 26 Mar 2026 20:04:48 -0400 Subject: [PATCH 08/12] ANALYSIS - IMPROVEMENT - Fix function name bug, add savefig/display to analyze scripts, increase margins and fix colorbar overlap Co-Authored-By: Claude Sonnet 4.6 --- .../analyze_example.jl | 10 +++- .../Solovev_ideal_example/analyze_example.jl | 10 +++- .../analyze_example.jl | 11 +++- src/Analysis/Equilibrium.jl | 12 ++-- src/Analysis/ForceFreeStates.jl | 22 +++---- src/Analysis/PerturbedEquilibrium.jl | 59 ++++++++++--------- 6 files changed, 70 insertions(+), 54 deletions(-) diff --git a/examples/DIIID-like_ideal_example/analyze_example.jl b/examples/DIIID-like_ideal_example/analyze_example.jl index dbf5da961..97d93bb87 100644 --- a/examples/DIIID-like_ideal_example/analyze_example.jl +++ b/examples/DIIID-like_ideal_example/analyze_example.jl @@ -1,12 +1,16 @@ using Pkg; Pkg.activate(joinpath(@__DIR__, "../..")) -using GeneralizedPerturbedEquilibrium, Plots, PlotlyJS +using GeneralizedPerturbedEquilibrium, Plots using GeneralizedPerturbedEquilibrium: Analysis -plotlyjs() +gr() h5path = joinpath(@__DIR__, "gpec.h5") # Summary plots p_eq = Analysis.Equilibrium.plot_equilibrium_summary(h5path) p_ffs = Analysis.ForceFreeStates.plot_ffs_summary(h5path) -p_pe = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path) \ No newline at end of file +p_pe = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path) + +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")) \ No newline at end of file diff --git a/examples/Solovev_ideal_example/analyze_example.jl b/examples/Solovev_ideal_example/analyze_example.jl index dbf5da961..97d93bb87 100644 --- a/examples/Solovev_ideal_example/analyze_example.jl +++ b/examples/Solovev_ideal_example/analyze_example.jl @@ -1,12 +1,16 @@ using Pkg; Pkg.activate(joinpath(@__DIR__, "../..")) -using GeneralizedPerturbedEquilibrium, Plots, PlotlyJS +using GeneralizedPerturbedEquilibrium, Plots using GeneralizedPerturbedEquilibrium: Analysis -plotlyjs() +gr() h5path = joinpath(@__DIR__, "gpec.h5") # Summary plots p_eq = Analysis.Equilibrium.plot_equilibrium_summary(h5path) p_ffs = Analysis.ForceFreeStates.plot_ffs_summary(h5path) -p_pe = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path) \ No newline at end of file +p_pe = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path) + +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")) \ No newline at end of file diff --git a/examples/Solovev_ideal_example_multi_n/analyze_example.jl b/examples/Solovev_ideal_example_multi_n/analyze_example.jl index ec1873cec..9e4ee2006 100644 --- a/examples/Solovev_ideal_example_multi_n/analyze_example.jl +++ b/examples/Solovev_ideal_example_multi_n/analyze_example.jl @@ -1,8 +1,8 @@ using Pkg; Pkg.activate(joinpath(@__DIR__, "../..")) -using GeneralizedPerturbedEquilibrium, Plots, PlotlyJS +using GeneralizedPerturbedEquilibrium, Plots using GeneralizedPerturbedEquilibrium: Analysis -plotlyjs() +gr() # Analyze the multi-n run h5path = joinpath(@__DIR__, "gpec.h5") @@ -11,8 +11,15 @@ p_eq = Analysis.Equilibrium.plot_equilibrium_summary(h5path) p_ffs = Analysis.ForceFreeStates.plot_ffs_summary(h5path) p_pe = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path) +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")) + # Analyze the single-n runs h5path_n1 = joinpath(@__DIR__, "single_n_1", "euler_n1.h5") h5path_n2 = joinpath(@__DIR__, "single_n_2", "euler_n2.h5") p_n1 = Analysis.ForceFreeStates.plot_ffs_summary(h5path_n1) p_n2 = Analysis.ForceFreeStates.plot_ffs_summary(h5path_n2) + +display(p_n1); Plots.savefig(p_n1, joinpath(@__DIR__, "single_n_1", "ffs_summary.png")) +display(p_n2); Plots.savefig(p_n2, joinpath(@__DIR__, "single_n_2", "ffs_summary.png")) diff --git a/src/Analysis/Equilibrium.jl b/src/Analysis/Equilibrium.jl index 66a2466e8..c82c419cb 100644 --- a/src/Analysis/Equilibrium.jl +++ b/src/Analysis/Equilibrium.jl @@ -41,7 +41,7 @@ function plot_qprofile(h5path; show_singular=true, save_path=nothing) title="", legend=false, xlims=(0, 1), - left_margin=5Plots.mm, + left_margin=10Plots.mm, bottom_margin=5Plots.mm, right_margin=8Plots.mm # annotations at x=1.0 need room on the right ) @@ -94,7 +94,7 @@ function plot_pressure_profile(h5path; save_path=nothing) title="", legend=false, xlims=(0, 1), - left_margin=5Plots.mm, + left_margin=10Plots.mm, bottom_margin=5Plots.mm ) for s in 1:msing @@ -135,7 +135,7 @@ function plot_f_profile(h5path; save_path=nothing) title="", legend=false, xlims=(0, 1), - left_margin=5Plots.mm, + left_margin=10Plots.mm, bottom_margin=5Plots.mm ) for s in 1:msing @@ -198,7 +198,7 @@ function plot_flux_surfaces(h5path; n_psi=11, n_theta=18, save_path=nothing) ylabel="Z [m]", aspect_ratio=:equal, legend=:outertopright, - left_margin=5Plots.mm, + left_margin=10Plots.mm, bottom_margin=5Plots.mm ) @@ -270,7 +270,7 @@ function plot_gse_by_theta(h5path; n_theta_lines=8, save_path=nothing) title="", yscale=:log10, xlims=(0, 1), - left_margin=5Plots.mm, + left_margin=10Plots.mm, bottom_margin=5Plots.mm ) for (i, itheta) in enumerate(theta_indices) @@ -326,7 +326,7 @@ function plot_gse_integrated(h5path; save_path=nothing) ylabel="log₁₀(integrated GSE)", title="Flux-surface-integrated Grad-Shafranov error", legend=false, - left_margin=5Plots.mm, + left_margin=10Plots.mm, bottom_margin=5Plots.mm ) diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index 610744437..8b806a682 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -46,7 +46,7 @@ function plot_mode_displacement(h5path; modes=1:5, save_path=nothing) xlabel="Norm. Poloidal Flux", ylabel="|ξ^ψ|", title=title_str, - left_margin=5Plots.mm, + left_margin=10Plots.mm, bottom_margin=5Plots.mm ) for m in modes @@ -60,7 +60,7 @@ end """ - plot_fixed_boundarystability_criterion(h5path; save_path=nothing) + plot_fixed_boundary_stability_criterion(h5path; save_path=nothing) Plot the stability criterion (smallest eigenvalue of W⁻¹, `crit`) vs ψ_N. A sign change in `crit` during integration indicates an ideal fixed-boundary instability. @@ -77,7 +77,7 @@ A sign change in `crit` during integration indicates an ideal fixed-boundary ins A `Plots.jl` plot object. """ -function plot_fixed_boundarystability_criterion(h5path; save_path=nothing) +function plot_fixed_boundary_stability_criterion(h5path; save_path=nothing) psi, crit = h5open(h5path, "r") do fid read(fid["integration/psi"]), read(fid["integration/crit"]) end @@ -89,7 +89,7 @@ function plot_fixed_boundarystability_criterion(h5path; save_path=nothing) ylabel="|Dᶜ|", title="Fixed-Boundary Stability", legend=false, - left_margin=5Plots.mm, + left_margin=10Plots.mm, bottom_margin=5Plots.mm ) hline!(p, [0.0]; linestyle=:dash, color=:black, label=nothing) @@ -144,8 +144,8 @@ function plot_energy_eigenvectors(h5path; matrix_type=:total, save_path=nothing) ylabel="Eigenmode Index", title="Total Energy Eigenvectors", colorbar_title="Harmonic Amplitude", - left_margin=5Plots.mm, - right_margin=10Plots.mm, + left_margin=10Plots.mm, + right_margin=20Plots.mm, bottom_margin=5Plots.mm ) @@ -194,7 +194,7 @@ function plot_eigenvalues(h5path; matrix_type=:total, save_path=nothing) legend=false, color=colors, markerstrokewidth=0, - left_margin=5Plots.mm, + left_margin=10Plots.mm, bottom_margin=5Plots.mm ) hline!(p, [0]; linestyle=:dash, color=:black, label=nothing) @@ -265,7 +265,7 @@ function plot_delta_prime(h5path; save_path=nothing) color=colors, markersize=7, markerstrokewidth=0, - left_margin=5Plots.mm, + left_margin=10Plots.mm, bottom_margin=5Plots.mm ) hline!(p, [0.0]; linestyle=:dash, color=:black, label=nothing) @@ -304,16 +304,16 @@ function plot_ffs_summary(h5path; save_path=nothing) haskey(fid, "vacuum/wt") && !isempty(read(fid["vacuum/wt"])) end - p_crit = plot_stability_criterion(h5path) + 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)) + 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)) + p = plot(p_crit, p_dp; layout=(1, 2), size=(1100, 500), right_margin=10Plots.mm) end isnothing(save_path) || savefig(p, save_path) diff --git a/src/Analysis/PerturbedEquilibrium.jl b/src/Analysis/PerturbedEquilibrium.jl index 6b2465a26..6ba24ddca 100644 --- a/src/Analysis/PerturbedEquilibrium.jl +++ b/src/Analysis/PerturbedEquilibrium.jl @@ -51,7 +51,7 @@ function plot_resonant_flux(h5path; save_path=nothing) p = plot(; xlabel="Norm. Poloidal Flux", ylabel="|Φ_res|", title="Resonant flux |Φ_res| per surface", legend=:outertopright, - left_margin=5Plots.mm, bottom_margin=5Plots.mm) + left_margin=10Plots.mm, bottom_margin=5Plots.mm) n_vals = unique(pe_n) for nn in n_vals @@ -109,7 +109,7 @@ function plot_island_widths(h5path; save_path=nothing) color=:steelblue, markersize=7, markerstrokewidth=0, - left_margin=5Plots.mm, + left_margin=10Plots.mm, bottom_margin=5Plots.mm ) for s in 1:msing @@ -164,7 +164,7 @@ function plot_chirikov_parameter(h5path; save_path=nothing) color=colors, markersize=7, markerstrokewidth=0, - left_margin=5Plots.mm, + left_margin=10Plots.mm, bottom_margin=5Plots.mm ) hline!(p, [1.0]; linestyle=:dash, color=:black, label=nothing) @@ -215,7 +215,7 @@ function plot_driven_delta_prime(h5path; save_path=nothing) p = plot(; xlabel="Norm. Poloidal Flux", ylabel="Re(Δ')", title="Tearing stability Δ' (PE)", legend=:outertopright, - left_margin=5Plots.mm, bottom_margin=5Plots.mm) + left_margin=10Plots.mm, bottom_margin=5Plots.mm) hline!(p, [0.0]; linestyle=:dash, color=:black, label=nothing) n_vals = unique(pe_n) @@ -301,7 +301,7 @@ function _plot_resonant_current(h5path) p = plot(; xlabel="Norm. Poloidal Flux", ylabel="|I_res|", title="Resonant current |I_res| per surface", legend=:outertopright, - left_margin=5Plots.mm, bottom_margin=5Plots.mm) + left_margin=10Plots.mm, bottom_margin=5Plots.mm) n_vals = unique(pe_n) for nn in n_vals @@ -382,7 +382,7 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) ylabel="|$(component)|", title="Mode spectrogram: $(component)", legend=:outertopright, - left_margin=5Plots.mm, + left_margin=10Plots.mm, bottom_margin=5Plots.mm ) cmap = cgrad(:roma, mpert; categorical=true) @@ -400,8 +400,8 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) ylabel="Norm. Poloidal Flux", title="", colorbar_title="|$(component)|", - left_margin=5Plots.mm, - right_margin=10Plots.mm, + left_margin=10Plots.mm, + right_margin=20Plots.mm, bottom_margin=5Plots.mm ) # Overlay rational surface locations as white dashed lines @@ -437,43 +437,44 @@ 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_island_widths(h5path) - p_energies = _plot_energies(h5path) + p_islands = plot_resonant_flux(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_energies, p_spectro; layout=l, size=(1100, 1100)) + p = plot(p_islands, p_bpsi, p_spectro; layout=l, size=(1100, 1100)) isnothing(save_path) || savefig(p, save_path) return p end -# Internal helper — energy breakdown bar chart -function _plot_energies(h5path) - key = "perturbed_equilibrium/energies/plasma_energy" - _has_pe_data(h5path, key) || - return plot(; title="No energy data", legend=false) +# Internal helper — |b_psi(m)| spectrum at the outermost psi surface +function _plot_bpsi_edge_spectrum(h5path) + base = "perturbed_equilibrium/response/" + _has_pe_data(h5path, base * "b_psi_real") || + return plot(; title="No b_psi data — run with perturbed equilibrium enabled", legend=false) - ep, ev, et = h5open(h5path, "r") do fid - read(fid["perturbed_equilibrium/energies/plasma_energy"]), - read(fid["perturbed_equilibrium/energies/vacuum_energy"]), - read(fid["perturbed_equilibrium/energies/total_energy"]) + data_r, data_i, mlow, mhigh = h5open(h5path, "r") do fid + read(fid[base * "b_psi_real"]), read(fid[base * "b_psi_imag"]), + read(fid["info/mlow"]), read(fid["info/mhigh"]) end - vals = [real(ep), real(ev), real(et)] # real() in case energies are stored as complex - names = ["Plasma", "Vacuum", "Total"] - colors = [:steelblue, :darkorange, :green] + mpert = mhigh - mlow + 1 + m_vals = mlow:mhigh + data = complex.(data_r, data_i) # (npsi, numpert_total) + bpsi_edge = abs.(data[end, 1:mpert]) p = bar( - names, vals; - ylabel="Energy", - title="Energy breakdown", + collect(m_vals), bpsi_edge; + xlabel="m", + ylabel="|b_ψ|", + title="b_ψ spectrum at ψ edge", legend=false, - color=colors, - left_margin=5Plots.mm, + color=:steelblue, + linewidth=0, + left_margin=10Plots.mm, bottom_margin=5Plots.mm ) - hline!(p, [0]; linestyle=:dash, color=:black, label=nothing) return p end From abc5aa8dc3729c0807a5bbe30e2513b6356884a3 Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 27 Mar 2026 09:45:37 -0400 Subject: [PATCH 09/12] ANALYSIS - FIX - Add PerturbedEquilibrium submodule to analysis docs Documenter was failing with missing_docs error because all PerturbedEquilibrium analysis functions had docstrings but no @autodocs block in analysis.md. Co-Authored-By: Claude Sonnet 4.6 --- docs/src/analysis.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/docs/src/analysis.md b/docs/src/analysis.md index 8940a3e4b..0e7f3c100 100644 --- a/docs/src/analysis.md +++ b/docs/src/analysis.md @@ -6,6 +6,7 @@ The Analysis module provides post-processing and visualization utilities for GPE - `ForceFreeStates`: Plotting functions for DCON-style ideal MHD stability results - `Equilibrium`: Plotting functions for equilibrium objects +- `PerturbedEquilibrium`: Plotting functions for perturbed equilibrium results ```@docs GeneralizedPerturbedEquilibrium.Analysis @@ -22,3 +23,9 @@ Modules = [GeneralizedPerturbedEquilibrium.Analysis.ForceFreeStates] ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.Analysis.Equilibrium] ``` + +## PerturbedEquilibrium + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Analysis.PerturbedEquilibrium] +``` From 20394986356c894a75fe7976852b8930fae65c5a Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 27 Mar 2026 09:47:33 -0400 Subject: [PATCH 10/12] DOCS - Add reminder to update docs coverage when adding new modules Co-Authored-By: Claude Sonnet 4.6 --- CLAUDE.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CLAUDE.md b/CLAUDE.md index d88b49ab3..c001177fb 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -478,6 +478,7 @@ This format is used for compiling release notes, so tags should be human-readabl - **Never remove packages from Project.toml** - If a package fails to load or resolve, run `Pkg.add(...)` or `Pkg.instantiate()` to fix the local environment. Do NOT remove the package from `Project.toml`. The developer works across multiple branches and machines, so environment drift is expected — the right fix is always to update the environment to satisfy the toml, not to trim the toml to match the current environment state. - **Indexing**: The codebase uses 0-based indexing in many places to match Fortran conventions, then converts to 1-based Julia indexing - **No step numbering in code comments** - Avoid annotations like "Step 1: do this" followed by "Step 2: do that". These get out of sync as code changes. Just describe the action without numbering. +- **Documentation coverage** - When adding a new module or submodule with public docstrings, add a corresponding `@autodocs` block in `docs/src/`. Documenter CI will fail with a `missing_docs` error if any exported docstring is not covered. The analysis submodule docs live in `docs/src/analysis.md`. - **Keep code comments concise** - A comment should be one line where possible. Do not write multi-line block comments explaining the current session's investigation, what was tried, what was wrong before, or why a specific file/path behaves differently. State what the code does and why at a general level. Example of too much detail: a 6-line block explaining that efit_by_inversion uses psilow>0 while CHEASE starts at 0, that the old code was removed, and that spline spikes result. Preferred: `# Replicate Fortran inverse.f: overwrite deta at axis (r²=0) by extrapolating from innermost surfaces.` ### Output Files From 571a9568b2e19bc0e2a5b7001a897d2b3db71924 Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 27 Mar 2026 09:56:16 -0400 Subject: [PATCH 11/12] ANALYSIS - IMPROVEMENT - Use plotlyjs() interactively, gr() for scripts Co-Authored-By: Claude Sonnet 4.6 --- examples/DIIID-like_ideal_example/analyze_example.jl | 2 +- examples/Solovev_ideal_example/analyze_example.jl | 2 +- examples/Solovev_ideal_example_multi_n/analyze_example.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/DIIID-like_ideal_example/analyze_example.jl b/examples/DIIID-like_ideal_example/analyze_example.jl index 97d93bb87..4c5a8eb7c 100644 --- a/examples/DIIID-like_ideal_example/analyze_example.jl +++ b/examples/DIIID-like_ideal_example/analyze_example.jl @@ -2,7 +2,7 @@ using Pkg; Pkg.activate(joinpath(@__DIR__, "../..")) using GeneralizedPerturbedEquilibrium, Plots using GeneralizedPerturbedEquilibrium: Analysis -gr() +isinteractive() ? plotlyjs() : gr() h5path = joinpath(@__DIR__, "gpec.h5") diff --git a/examples/Solovev_ideal_example/analyze_example.jl b/examples/Solovev_ideal_example/analyze_example.jl index 97d93bb87..4c5a8eb7c 100644 --- a/examples/Solovev_ideal_example/analyze_example.jl +++ b/examples/Solovev_ideal_example/analyze_example.jl @@ -2,7 +2,7 @@ using Pkg; Pkg.activate(joinpath(@__DIR__, "../..")) using GeneralizedPerturbedEquilibrium, Plots using GeneralizedPerturbedEquilibrium: Analysis -gr() +isinteractive() ? plotlyjs() : gr() h5path = joinpath(@__DIR__, "gpec.h5") diff --git a/examples/Solovev_ideal_example_multi_n/analyze_example.jl b/examples/Solovev_ideal_example_multi_n/analyze_example.jl index 9e4ee2006..3a1e46134 100644 --- a/examples/Solovev_ideal_example_multi_n/analyze_example.jl +++ b/examples/Solovev_ideal_example_multi_n/analyze_example.jl @@ -2,7 +2,7 @@ using Pkg; Pkg.activate(joinpath(@__DIR__, "../..")) using GeneralizedPerturbedEquilibrium, Plots using GeneralizedPerturbedEquilibrium: Analysis -gr() +isinteractive() ? plotlyjs() : gr() # Analyze the multi-n run h5path = joinpath(@__DIR__, "gpec.h5") From 0adbc17f2d57ad8d4dfb199c325783de5797e854 Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 27 Mar 2026 09:56:51 -0400 Subject: [PATCH 12/12] Examples - MINOR - Returns all to save_interval 3 --- examples/DIIID-like_ideal_example/gpec.toml | 2 +- examples/Solovev_ideal_example_3D/gpec.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 25c660829..12f073263 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -60,7 +60,7 @@ ion_flag = true # Include ion dW_k when kin_flag is true electron_flag = false # Include electron dW_k when kin_flag is true eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations -save_interval = 1 # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. +save_interval = 3 # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e4 # Maximum fraction of solutions allowed before re-normalized diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml index d9e8b6635..01961b4bc 100644 --- a/examples/Solovev_ideal_example_3D/gpec.toml +++ b/examples/Solovev_ideal_example_3D/gpec.toml @@ -52,7 +52,7 @@ eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e3 # Maximum fraction of solutions allowed before re-normalized force_wv_symmetry = true # Forces vacuum energy matrix symmetry -save_interval = 10 # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. +save_interval = 3 # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. [Wall] shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath)