ANALYSIS - Add equilibrium, stability, and PE visualization toolkit

CLAUDE.md

@@ -475,8 +475,10 @@ 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.
+- **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

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"

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]
+```

examples/DIIID-like_ideal_example/analyze_example.jl

@@ -1,10 +1,16 @@
 using Pkg;
 Pkg.activate(joinpath(@__DIR__, "../.."))
-using GeneralizedPerturbedEquilibrium, Plots, PlotlyJS
+using GeneralizedPerturbedEquilibrium, Plots
 using GeneralizedPerturbedEquilibrium: Analysis
-plotlyjs()
+isinteractive() ? plotlyjs() : gr()
 
 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_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"))

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 = 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
 
 [ForcingTerms]
 forcing_data_file = "forcing.dat"       # Path to forcing data file (n, m, complex amplitude)

examples/Solovev_ideal_example/analyze_example.jl

@@ -1,10 +1,16 @@
 using Pkg;
 Pkg.activate(joinpath(@__DIR__, "../.."))
-using GeneralizedPerturbedEquilibrium, Plots, PlotlyJS
+using GeneralizedPerturbedEquilibrium, Plots
 using GeneralizedPerturbedEquilibrium: Analysis
-plotlyjs()
+isinteractive() ? plotlyjs() : gr()
 
 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_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"))

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)

examples/Solovev_ideal_example_multi_n/analyze_example.jl

@@ -1,17 +1,25 @@
 using Pkg;
 Pkg.activate(joinpath(@__DIR__, "../.."))
-using GeneralizedPerturbedEquilibrium, Plots, PlotlyJS
+using GeneralizedPerturbedEquilibrium, Plots
 using GeneralizedPerturbedEquilibrium: Analysis
-plotlyjs()
+isinteractive() ? plotlyjs() : gr()
 
 # 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_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_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)
+
+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"))

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