This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
GPEC (Generalized Perturbed Equilibrium Code, Julia implementation) is a comprehensive Julia reimplementation of the GPEC suite for MHD analysis of fusion plasmas. The code performs equilibrium reconstruction, ideal MHD stability analysis, and perturbed equilibrium calculations including plasma response and singular surface coupling diagnostics.
Relationship to Fortran GPEC: This Julia GPEC is an evolution of the Fortran GPEC code suite, available at https://github.com/PrincetonUniversity/GPEC. When users reference "the Fortran code", "the original GPEC", or "Fortran GPEC", they are referring to that Fortran codebase. This Julia implementation reimplements and extends GPEC's functionality with improved performance and maintainability.
Local GPEC Repository: For code conversion or comparison with the original Fortran implementation, check for a local GPEC repository at ~/Code/gpec. If not found at this location, ask the user for the correct path.
This codebase is implemented in Julia. References to “Fortran GPEC” or “legacy VACUUM” refer to the original upstream Fortran codebase, not a runtime dependency of this Julia implementation.
Current Development Focus: The perturbed_equilibrium branch is implementing full GPEC-style perturbed equilibrium functionality, including singular coupling analysis, island formation diagnostics, and mode-space field reconstruction.
IMPORTANT: The papers in docs/resources/ provide the theoretical foundation for GPEC's algorithms and should be referenced to understand what the code is doing. Citing equations from these papers in code comments and annotations is strongly encouraged to maintain traceability between theory and implementation.
The Vacuum module implements the methods described in:
-
Chance et al. (1997): "Vacuum calculations in azimuthally symmetric geometry"
- Location:
docs/resources/1997-Chance-Vacuum_calculations_in_azimuthally_symmetric_geometry.pdf - Published: Physics of Plasmas 4, 2161 (1997)
- Link: https://pubs.aip.org/aip/pop/article-abstract/4/6/2161/263192
- Describes: Fundamental vacuum response calculation method for tokamak geometry
- Location:
-
Chance et al. (2007): "Calculation of the vacuum Green's function valid even for high toroidal mode numbers in tokamaks"
- Location:
docs/resources/2007-Chance-Calculation of the vacuum Greens function valid even for high toroidal mode numbers in tokamaks.pdf - Published: Physics of Plasmas 14, 052506 (2007)
- Describes: Improved Green's function calculation for high-n modes
- Location:
The ForceFreeStates module (ideal MHD stability analysis) implements methods from:
-
Glasser (2016): "The direct criterion of Newcomb for the ideal MHD stability of an axisymmetric toroidal plasma"
- Location:
docs/resources/2016-Glasser-The_direct_criterion_of_Newcomb_for_the_ideal_MHD_stability_of_an_axisymmetric_toroidal_plasma.pdf - Published: Physics of Plasmas 23, 112506 (2016)
- Describes: FUNDAMENTAL PAPER - Newcomb's criterion for ideal MHD stability (current implementation)
- Location:
-
Glasser (2018): "A Riccati solution for the ideal MHD plasma response with applications to real-time stability control"
- Location:
docs/resources/2018-Glasser-A Riccati solution for the ideal MHD plasma response with applications to real-time stability control.pdf - Published: Physics of Plasmas 25, 032507 (2018)
- Describes: Riccati method for ideal MHD eigenvalue problem
- Location:
The PerturbedEquilibrium module implements GPEC-style perturbed equilibrium calculations from:
-
Park et al. (2007a): "Computation of three-dimensional tokamak and spherical torus equilibria"
- Location:
docs/resources/2007-Park-Computation_of_three-dimensional_tokamak_and_spherical_torus_equilibria-compressed.pdf - Published: Physics of Plasmas 14, 052110 (2007)
- Describes: 3D equilibrium perturbations in toroidal geometry
- Location:
-
Park et al. (2007b): "Control of Asymmetric Magnetic Perturbations in Tokamaks"
- Location:
docs/resources/2007-Park-Control_of_Asymmetric_Magnetic_Perturbations_in_Tokamaks.pdf - Published: Physical Review Letters 99, 195003 (2007)
- Describes: Plasma response to resonant magnetic perturbations (RMP)
- Location:
-
Park et al. (2009): "Importance of plasma response to nonaxisymmetric perturbations in tokamaks"
- Location:
docs/resources/2009-Park-Importance_of_plasma_response_to_nonaxisymmetric_perturbations_in_tokamaks-compressed.pdf - Published: Physics of Plasmas 16, 056115 (2009)
- Describes: Self-consistent plasma response calculation
- Location:
-
Park et al. (2011): "Kinetic energy principle and neoclassical toroidal torque in tokamaks"
- Location:
docs/resources/2011-Park-Physics_of_Plasmas_Kinetic_energy_principle_and_neoclassical_toroidal_torque_in_tokamaks.pdf - Published: Physics of Plasmas 18, 110702 (2011)
- Describes: Energy principle for perturbed equilibria
- Location:
-
Park et al. (2017): "Self-consistent perturbed equilibrium with neoclassical toroidal torque in tokamaks"
- Location:
docs/resources/2017-Park-Self_consistent_perturbed_equilibrium_with_neoclassical_toroidal_torque_in_toka.pdf - Published: Physics of Plasmas 24, 032505 (2017)
- Describes: Self-consistent coupling with neoclassical effects
- Location:
GPEC will eventually implement resistive MHD stability analysis based on:
-
Glasser (2016): "Computation of resistive instabilities by matched asymptotic expansions"
- Location:
docs/resources/2016-Glasser-Computation_of_resistive_instabilities_by_matched_asymptotic_expansions-compressed.pdf - Published: Physics of Plasmas 23, 072505 (2016)
- Describes: Resistive stability analysis and Δ' calculation via matched asymptotic expansions
- Location:
-
Glasser (2018): "A robust solution for the resistive MHD toroidal Δ′ matrix in near real-time"
- Location:
docs/resources/2018-Glasser-A robust solution for the resistive MHD toroidal Delta-prime matrix in near real-time.pdf - Published: Physics of Plasmas 25, 032501 (2018)
- Describes: Fast computation of Δ' matrix for resistive stability
- Location:
-
Wang et al. (2020): "Modeling of resistive plasma response in toroidal geometry using an asymptotic matching approach"
- Location:
docs/resources/2020-Wang-Modeling of resistive plasma response in toroidal geometry using an asymptotic matching approach.pdf - Published: Physics of Plasmas 27, 122509 (2020)
- Describes: Asymptotic matching for resistive plasma response
- Location:
GPEC will eventually port the PENTRC (Perturbed Equilibrium Neoclassical Toroidal viscosity in Realistic geometry Code) functionality from the Fortran GPEC suite. This is described in:
-
Logan & Park (2013): "Neoclassical toroidal viscosity in perturbed equilibria with general tokamak geometry"
- Location:
docs/resources/2013-Logan-Neoclassical_toroidal_viscosity_in_perturbed_equilibria_with_general_tokamak_geometry.pdf - Published: Physics of Plasmas 20, 122507 (2013)
- Describes: Neoclassical toroidal viscosity (NTV) in perturbed equilibria
- Location:
-
Logan (2015): "Electromagnetic Torque in Tokamaks with Toroidal Asymmetries"
- Location:
docs/resources/2015-Logan-Electromagnetic_Torque_in_Tokamaks_with_Toroidal_Asymmetries-compressed.pdf - Published: PhD Thesis, Princeton University (2015)
- Describes: Complete PENTRC theory and implementation
- Location:
- Various (2009): "Nonambipolar Transport by Trapped Particles in Tokamaks"
- Location:
docs/resources/2009-Nonambipolar_Transport_by_Trapped_Particles_in_Tokamaks.pdf - Describes: Neoclassical transport theory
- Location:
# Run all tests
julia --project=. -e 'using Pkg; Pkg.activate("."); Pkg.instantiate(); include("test/runtests.jl")'
# Run specific test file
julia --project=. test/runtests.jl test/runtests_solovev.jl
# Available test files:
# - test/runtests_vacuum_julia.jl # Julia vacuum module
# - test/runtests_solovev.jl # Analytical equilibrium
# - test/runtests_ode.jl # ODE integration
# - test/runtests_sing.jl # Singular surface handling
# - test/runtests_fullruns.jl # End-to-end tests# Build documentation locally
julia --project=. build_docs_local.jl
# Documentation hosted at: https://openfusiontoolkit.github.io/GPEC/dev/For faster recompilation during development, use Revise.jl (installed in global environment, not in Project.toml):
using Revise
using GeneralizedPerturbedEquilibriumUse the generic benchmarking tool at benchmarks/benchmark_git_branches.jl to compare performance between branches or commits.
Tool usage:
# Compare feature branch against develop
julia benchmarks/benchmark_git_branches.jl \
--example examples/DIIID-like_ideal_example \
--branch1 develop \
--branch2 feature-branch
# Compare specific commits
julia benchmarks/benchmark_git_branches.jl \
--example examples/DIIID-like_ideal_example \
--commit1 abc123 \
--commit2 def456
# Compare current develop vs develop from 1 month ago
julia benchmarks/benchmark_git_branches.jl \
--example examples/DIIID-like_ideal_example \
--branch1 develop \
--commit1 HEAD~10 \
--branch2 developDefault benchmark case: examples/DIIID-like_ideal_example
Reported metrics:
- Eigenmode energy (
et[1]) - First eigenvalue; verifies calculation correctness - Integration steps - Total ODE solver steps
- Runtime (warmed) - Wall-clock time averaged over multiple warm runs (JIT warmup handled automatically)
- Commit hash - Git commit of code tested
The tool automatically:
- Handles JIT warmup (runs example 3 times, averages last 2)
- Switches between branches/commits
- Stashes uncommitted changes if necessary
- Restores original branch when done
- Reports comparison with percentage differences
Important notes:
- Working directory should be clean or changes will be stashed during branch switching
- Tool requires HDF5.jl for reading euler.h5 output
- Each benchmark run takes several minutes per branch (includes compilation + warm runs)
GPEC follows a three-stage analysis pipeline:
- Equilibrium → Solve Grad-Shafranov equation, compute flux surfaces, safety factor q-profile
- Stability Analysis → Solve ideal MHD eigenvalue problem (DCON-style), identify singular surfaces
- Perturbed Equilibrium → Compute plasma response to external fields, analyze singular coupling and island formation
This workflow is reflected in the modular structure and data flow.
GPEC consists of seven main modules organized in src/:
-
Splines (
src/Splines/) - Numerical interpolation libraryCubicSpline.jl- 1D cubic spline interpolationBicubicSpline.jl- 2D bicubic spline interpolationFourierSpline.jl- Fourier-based spline interpolation- Status: Mature, pure Julia implementation
-
Utilities (
src/Utilities/) - Shared computational toolsFourierTransforms.jl- Efficient Fourier transform utilities with pre-computed basis functions- Provides type-stable functor pattern for repeated transforms
- Used by Vacuum and PerturbedEquilibrium modules
-
Equilibrium (
src/Equilibrium/) - MHD equilibrium solvers- Main entry point:
setup_equilibrium(path)orsetup_equilibrium(config) - Supports multiple equilibrium types:
efit- EFIT g-file formatchease,chease2- CHEASE equilibrium code formatslar- Large Aspect Ratio analytical modelsol- Solovev analytical equilibrium
- Key files:
EquilibriumTypes.jl- Core data structuresReadEquilibrium.jl- Parsing equilibrium filesDirectEquilibrium.jl- Direct Grad-Shafranov solverInverseEquilibrium.jl- Inverse equilibrium solverAnalyticEquilibrium.jl- Analytical solutions
- Status: Stable and feature-complete
- Main entry point:
-
Vacuum (
src/Vacuum/) - Vacuum field calculations and Green's functions- Computes vacuum response matrices for ideal MHD analysis
- Calculates both interior (grri) and exterior (grre) Green's functions
- Main functions:
compute_vacuum_response()- Pure Julia implementation
- Key files:
VacuumStructs.jl- Data structuresVacuumInternals.jl- Core algorithmsVacuumFromEquilibrium.jl- Integration with equilibrium data
- Status: Pure Julia implementation complete and available
-
ForceFreeStates (
src/ForceFreeStates/) - Ideal MHD stability analysis (DCON-style)- Solves ideal MHD eigenvalue problem with force-free boundary conditions
- Identifies singular surfaces where ξ·∇ψ = 0
- Key files:
ForceFreeStatesStructs.jl- Core data structuresOde.jl- ODE solver for Euler-Lagrange equationsSing.jl- Singular point handling and layer analysisFourfit.jl- Fourier fitting routinesFixedBoundaryStability.jl- Fixed boundary analysisFree.jl- Free boundary stabilityMercier.jl- Mercier stability criterion
- Status: Stable, core DCON functionality implemented
-
ForcingTerms (
src/ForcingTerms/) - External field specification- Handles external magnetic field perturbations (coils, RMP, etc.)
- Supports ASCII and HDF5 forcing data formats
ForcingModedata structure specifies amplitude and phase for each (m,n) component- Status: Complete and functional
-
PerturbedEquilibrium (
src/PerturbedEquilibrium/) - GPEC-style plasma response- Computes plasma response to external forcing
- Calculates singular coupling metrics at rational surfaces
- Key files:
PerturbedEquilibrium.jl- Main entry pointPerturbedEquilibriumStructs.jl- Data structuresResponseMatrices.jl- Permeability matrix calculationFieldReconstruction.jl- Mode-space field reconstructionResponse.jl- Plasma response computationSingularCoupling.jl- Singular surface analysis including:- Delta prime (Δ') tearing stability parameter
- Resonant flux and currents at rational surfaces
- Island half-widths and Chirikov parameters
- Green's functions at interior flux surfaces
- Surface inductance for singular surfaces
Utils.jl- Helper functions
- Status: Active development on
perturbed_equilibriumbranch
Unified Configuration File: gpec.toml
All GPEC modules are configured via a single TOML file with the following sections:
[Equilibrium]- Equilibrium solver settings[Wall]- Wall geometry and vacuum region[ForceFreeStates]- Stability analysis parameters[PerturbedEquilibrium]- Perturbed equilibrium settings[ForcingTerms]- External field specification
Key parameters:
force_termination- Set totrueto exit after equilibrium/stability (skip perturbed equilibrium)output_file- Output filename (default:gpec.h5)
Example configuration files are provided in:
examples/Solovev_ideal_example/gpec.tomlexamples/DIIID-like_ideal_example/gpec.toml
Note: Legacy configuration files (equil.toml, vac.in) are deprecated.
The complete GPEC analysis pipeline:
-
Equilibrium Setup:
setup_equilibrium(config)reads configuration fromgpec.toml- Parses equilibrium data (EFIT, CHEASE, or analytical)
- Runs Grad-Shafranov solver (direct or inverse)
- Computes global parameters: q-profile, pressure, current density, β
- Creates bicubic splines for (ψ, θ, φ) → (R, Z, Φ) mapping
- Outputs:
PlasmaEquilibriumobject
-
Vacuum Response:
- Initialize plasma and wall surfaces from equilibrium
- Compute vacuum response matrices (wv, grri, grre)
- Calculate both interior and exterior Green's functions
- Pure Julia implementation
-
Stability Analysis (ForceFreeStates):
- Solve ideal MHD Euler-Lagrange equations via ODE integration
- Identify singular surfaces where q = m/n
- Compute Δ' at each singular surface
- Calculate potential and kinetic energies
- Check Mercier and ballooning stability criteria
- Outputs: Eigenmode structure ξ(ψ,θ)
-
Perturbed Equilibrium (GPEC-style):
- Load external forcing data (coil fields, RMP configuration)
- Compute plasma response using permeability matrices
- Reconstruct mode-space fields (ξ_modes, b_modes)
- Calculate singular coupling metrics at rational surfaces:
- Δ' (tearing stability parameter)
- Island half-widths
- Chirikov overlap parameter
- Resonant flux and currents
- Outputs:
PerturbedEquilibriumStatewith response fields and diagnostics
-
Output:
- All results saved to single HDF5 file (default:
gpec.h5) - HDF5 groups:
input/,info/,equil/,splines/,locstab/,integration/,singular/,vacuum/, and perturbed equilibrium data
- All results saved to single HDF5 file (default:
PlasmaEquilibrium- Main equilibrium container with bicubic splines (rzphi), 1D profiles (sq), and global parametersEquilibriumConfig- Configuration loaded from TOML files
VacuumInput- Input parameters for vacuum calculationsWallShapeSettings- Wall geometry configuration
SingType- Singular surface data including:- Rational surface location (ψ, q = m/n)
- Δ' (tearing stability parameter)
- Eigenmode structure at singular surface
- NEW: Green's functions (grri, grre) at interior singular surfaces
- Surface inductance
PerturbedEquilibriumControl- User-facing TOML configuration parametersPerturbedEquilibriumInternal- Internal state with mode arraysPerturbedEquilibriumState- Results including:- Response fields (ξ_modes, b_modes) in mode space
- Singular coupling matrices [msing × numpert_total]
- Island diagnostics (half-widths, Chirikov parameters)
ForcingMode- External forcing specification (m, n, amplitude, phase)
GeneralizedPerturbedEquilibrium
├── Splines (foundation)
├── Utilities (shared tools)
│ └── FourierTransforms
├── Equilibrium (uses Splines)
├── Vacuum (uses Splines, Equilibrium, Utilities)
├── ForcingTerms (data I/O)
├── ForceFreeStates (uses Equilibrium, Vacuum, Splines)
└── PerturbedEquilibrium (uses ForceFreeStates, Vacuum, ForcingTerms, Utilities)
This project uses GitFlow:
- Two permanent branches:
mainanddevelop mainbranch updated only at release-ready stagesdevelopbranch for integration of features- Feature branches off
develop, merged back with--no-ff
Current Development:
- Active branch:
perturbed_equilibrium- Major feature implementing GPEC-style perturbed equilibrium calculations - Previous work:
vacuum_juliabranch (merged) - Pure Julia vacuum implementation
CODE - TAG - Detailed message
Where:
- CODE: Module name (EQUIL, VAC, VACUUM, ForceFreeStates, PERTURBED EQUILIBRIUM, etc.)
- TAG: Type descriptor (WIP, MINOR, IMPROVEMENT, BUG FIX, NEW FEATURE, REFACTOR, CLEANUP, etc.)
Examples:
PERTURBED EQUILIBRIUM - NEW FEATURE - Implement singular coupling diagnosticsVAC - IMPROVEMENT - Add dual Green's function computationEQUIL - BUG FIX - Fixed separatrix finding for high kappaForceFreeStates - REFACTOR - Unified singular surface data structure
This format is used for compiling release notes, so tags should be human-readable and descriptive.
- Julia version: 1.11 is the target version
- 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.
- Default output:
gpec.h5(previouslyeuler.h5in older versions) - Legacy diagnostic files: When modifying equilibrium code, remember that older versions output
gsec.h5,gse.h5,gsei.h5- these are now consolidated intogpec.h5
- Perturbed equilibrium module: Active development of GPEC-style singular coupling analysis
- Configuration: All settings now in unified
gpec.tomlfile
- Pure Julia implementations are available for all major components and offer comparable or better performance than Fortran
- Benchmarks available in
benchmark/directory for Fourier transforms and vacuum calculations - Pre-commit hooks are configured for notebook cleaning and Julia formatting (see
docs/src/set_up.mdfor developer setup)
- When resolving git conflicts, do not simply accept one side.
- Analyze what each side changed and WHY before producing a resolution.
- Produce a merged version incorporating both sets of changes.
- If both sides renamed the same symbol differently, prefer the current (ours) branch convention.
- When a rename on one side conflicts with a logic change on the other, apply the logic change using the renamed symbol.
- If a conflict involves changes to numerical parameters (tolerances, boundary conditions, grid sizes), flag for human review rather than guessing.
- Flag any conflicts where the combination is ambiguous for human review.