Equilibrium Module
The Equilibrium module provides tools to read, construct, and analyze magnetohydrodynamic (MHD) plasma equilibria. It wraps a variety of equilibrium file formats (EFIT, CHEASE, SOL, LAR, and others), runs the appropriate direct or inverse solver, and returns a processed PlasmaEquilibrium object ready for downstream calculations.
Overview
Key responsibilities of the module:
- Read equilibrium input files and TOML configuration (see
EquilibriumConfig). - Provide convenient constructors for analytic / model equilibria (Large
Aspect Ratio, Solovev).- Build spline representations used throughout the code (1D cubic and
2D bicubic splines).- Run the direct or inverse equilibrium solver and post-process results
(global parameters, q-profile finding, separatrix finding, GSE checks).The module exposes a small public API that covers setup, configuration, and common analyses used by other GPEC components (e.g. ForceFreeStates, vacuum interfaces).
EFIT solver strategies
When eq_type is one of the EFIT-based options, three solver strategies are available:
eq_type | Method | Best for |
|---|---|---|
"efit" | Geometric-angle field-line ODE | Default; fast and robust for standard cases |
"efit_arclength" | Arc-length field-line ODE | Near-separatrix surfaces where the geometric-angle ODE becomes singular |
"efit_by_inversion" | Contour.jl marching-squares → inverse solver | Highest geometric accuracy; avoids ODE singularities entirely |
efit integrates field lines using the geometric angle η (polar angle from the magnetic axis, 0→2π) as the independent variable. Position at each step is computed as R = R₀ + rfac·cos(η), so the RHS denominator is Bz·cos(η) − Br·sin(η) — the projection of B_pol onto the radial direction. Near an X-point where Bp → 0 this denominator passes through zero, causing a coordinate singularity: the ODE steps become extremely small or the solver fails to converge.
efit_arclength uses arc length along the flux surface as the independent variable instead. The tangent direction (dR/ds, dZ/ds) = ∇ψ⊥/|∇ψ| has unit magnitude by construction, so the position equations have no denominator and remain well-behaved all the way to the separatrix. The 1/Bp factors in the accumulated integrals still diverge near X-points, but their solver tolerances are set large so they do not restrict the step size — only the position tracking drives adaptivity. The two methods produce identical results when both succeed; efit_arclength extends the reliable range of psihigh closer to 1.
efit_by_inversion traces flux surface level sets directly from ψ(R,Z) using marching squares, resamples each closed curve to a uniform geometric-angle grid, and feeds the result into the inverse equilibrium solver (the same path used by CHEASE input). The Cartesian evaluation grid is clipped to the separatrix bounding box and its resolution is set adaptively from a bilinear interpolation error bound, so no manual tuning is needed.
Radial grid packing
The default grid_type = "log_asymptotic" uses a three-region grid that respects the asymptotic behavior of q near both the magnetic axis and the separatrix:
- Core (ψ < 0.15): geometric spacing in log(ψ) — handles the axis where profiles behave as ψⁿ
- Middle (0.15 ≤ ψ ≤ 0.95): uniform spacing — preserves resolution through the pedestal region
- Edge (ψ > 0.95): geometric spacing in log(1−ψ) — tracks the logarithmic divergence q ~ −A·ln(1−ψ) near a diverted separatrix
With mpsi = 0 (the default), the number of radial knots is chosen automatically from the psi_accuracy parameter (target absolute error in q). Two probe field-line integrations near psihigh estimate the local log-slope A, and the knot count is set so the cubic spline error stays below psi_accuracy throughout the domain. The legacy grid_type = "ldp" (sin²-spaced) and explicit mpsi are still supported.
API Reference
GeneralizedPerturbedEquilibrium.Equilibrium.AnalyticEqSpec — Type
AnalyticEqSpecSingle-source-of-truth entry describing one analytic-equilibrium kind: the gpec.toml section that carries its parameters, the *Config type built from that section, and the run function that turns the config into solver input. tj_analytic and tj_analytic_direct share TJAnalyticConfig but bind different run functions, so run_fn is per-entry rather than derived from config_type.
Fields
section::String—gpec.tomlsection key (e.g."SOL_INPUT")config_type::DataType—*Configtype, constructed via its(::Dict)/(::String)ctorrun_fn::Function—(EquilibriumConfig, config) -> eq_inputsolver entry point
GeneralizedPerturbedEquilibrium.Equilibrium.DirectBField — Type
DirectBFieldInternal mutable struct to hold B-field components and their derivatives at a point. It is used as a temporary workspace to avoid allocations in tight loops.
GeneralizedPerturbedEquilibrium.Equilibrium.DirectIngest — Type
DirectIngestThe serializable raw arrays and scalars captured by a direct-equilibrium reader (read_efit, read_imas, sol_run is analytic and does not use this) — everything needed to rebuild a DirectRunInput's splines without re-reading the original g-file. Stored on DirectRunInput.ingest / PlasmaEquilibrium.ingest and dumped to gpec.h5 so a run can be replayed from any directory. The interpolants themselves are not serializable; they are reconstructed from these nodes by build_direct_from_ingest.
Fields
sq_xs::Vector{Float64}— normalized-ψ knots for the 1D profile splinesq_fs::Matrix{Float64}— 1D profile node values (F, μ₀P, q, √ψ_norm)psi_xs::Vector{Float64}— R grid for the 2D flux map [m]psi_ys::Vector{Float64}— Z grid for the 2D flux map [m]psi_rz::Matrix{Float64}— processed poloidal flux on the (R, Z) grid [Wb/rad]rmin/rmax/zmin/zmax::Float64— computational-grid bounds [m]psio::Float64— total flux difference |ψaxis - ψboundary| [Wb/rad]bt_sign::Int— sign of the toroidal field (+1 or -1)
GeneralizedPerturbedEquilibrium.Equilibrium.DirectRunInput — Type
DirectRunInput(...)A container struct that bundles all necessary inputs for the direct_run function. It is created by one of the equilibrium file read-in functions after processing the raw equilibrium data and preparing the initial splines.
Fields
config::EquilibriumConfigThe equilibrium configuration object.sq_in1D spline data versus normalized poloidal fluxpsin. Quantities:F = R * B_t— toroidal flux function [m·T]μ₀ * Pressure— plasma pressure (non-negative) [T²]q— safety factor profile√ψ_norm— square root of normalized flux
psi_in2D cubic interpolant on the (R, Z) grid [m]. The values correspond to the poloidal flux adjusted to be zero at the boundary [Wb/rad]. Definitions:ψ(R, Z) = ψ_boundary - ψ(R, Z)ψ = ψ * sign(ψ(centerR, centerZ))- 1D profiles are represented by
CubicInterpolantorCubicSeriesInterpolant - 2D flux surfaces by
CubicInterpolantND
- 1D profiles are represented by
psi_in_xs::Vector{Float64}— R coordinate grid for psi_in [m]psi_in_ys::Vector{Float64}— Z coordinate grid for psi_in [m]rmin::Float64— Minimum R-coordinate of the computational grid [m]rmax::Float64— Maximum R-coordinate of the computational grid [m]zmin::Float64— Minimum Z-coordinate of the computational grid [m]zmax::Float64— Maximum Z-coordinate of the computational grid [m]psio::Float64— Total flux difference|ψ_axis - ψ_boundary|[Wb/rad]bt_sign::Int— Sign of the toroidal field (+1 or -1); read from fpol sign in EFIT g-filesingest::EquilibriumIngest— captured raw arrays for thegpec.h5rerun snapshot (aDirectIngestfor file-based reads, ornothingfor analytic equilibria)
GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig — Type
Outer constructor for EquilibriumConfig from a parsed TOML dictionary
GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig — Type
EquilibriumConfig(...)A mutable struct containing configuration parameters for equilibrium reconstruction. Bundles all necessary settings originally specified in the equil fortran namelists.
Fields
eq_type::String- Type of equilibrium file ("efit", "solovev", "lar", etc.)eq_filename::String- Path to equilibrium input filejac_type::String- Jacobian coordinate type ("hamada", "pest", "equal_arc", "boozer", "park", "other")power_bp::Int- Poloidal field power exponent for Jacobianpower_b::Int- Total field power exponent for Jacobianpower_r::Int- Major radius power exponent for Jacobianpower_rc::Int- Minor radius (rfac = √((R-R₀)²+(Z-Z₀)²)) power exponent for Jacobianr0exp::Float64- Major radius normalization for CHEASE/EQDSK [m]b0exp::Float64- On-axis toroidal field normalization for CHEASE/EQDSK [T]grid_type::String- Grid type for flux surface discretization ("log_asymptotic", "ldp", "pow1")psilow::Float64- Lower limit of normalized flux coordinatepsihigh::Float64- Upper limit of normalized flux coordinatempsi::Int- Number of radial grid points (0 = auto-compute from psi_accuracy)psi_accuracy::Float64- Target absolute error in q for auto-mpsi (used when mpsi=0 and gridtype="logasymptotic")mtheta::Int- Number of poloidal grid pointsnewq0::Int- Override for on-axis safety factor (0 = use input value)etol::Float64- Error tolerance for equilibrium solverforce_termination::Bool- Terminate after equilibrium setup (skip stability calculations)use_galgrid::Bool- Use the same grid as galerkin method
GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig — Method
Outer constructor for EquilibriumConfig that enables a toml file interface for specifying the configuration settings
DEPRECATED: Use [Equilibrium] section in gpec.toml instead
GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumParameters — Type
EquilibriumParametersA mutable struct containing computed equilibrium parameters and diagnostic flags.
Fields
ro::Union{Nothing,Float64}- R-coordinate of the magnetic axis [m]zo::Union{Nothing,Float64}- Z-coordinate of the magnetic axis [m]psio::Union{Nothing,Float64}- Total flux difference |ψaxis - ψboundary| [Wb/rad]rsep::Union{Nothing,Vector{Float64}}- R-coordinates of the plasma boundary [m]zsep::Union{Nothing,Vector{Float64}}- Z-coordinates of the plasma boundary [m]rext::Union{Nothing,Vector{Float64}}- R-coordinates of the plasma edge [m]zext::Union{Nothing,Vector{Float64}}- Z-coordinates of the plasma edge [m]psi0::Union{Nothing,Float64}- Normalized poloidal flux at reference locationb0::Union{Nothing,Float64}- Total magnetic field strength at the axis [T]q0::Union{Nothing,Float64}- Safety factor at the axisqmin::Union{Nothing,Float64}- Minimum safety factor in the plasmaqmax::Union{Nothing,Float64}- Maximum safety factor in the plasmaqa::Union{Nothing,Float64}- Safety factor at the plasma edgeq95::Union{Nothing,Float64}- Safety factor at 95% flux surfaceqextrema_psi::Union{Nothing,Vector{Float64}}- Normalized flux values at q extremaqextrema_q::Union{Nothing,Vector{Float64}}- Safety factor values at extremamextrema::Union{Nothing,Int}- Number of extrema in q-profilepsi_norm::Union{Nothing,Float64}- Normalized poloidal fluxb_norm::Union{Nothing,Float64}- Normalized magnetic field strengthpsi_axis::Union{Nothing,Float64}- Poloidal flux at the axispsi_boundary::Union{Nothing,Float64}- Poloidal flux at the boundarypsi_boundary_norm::Union{Nothing,Float64}- Normalized boundary fluxpsi_axis_norm::Union{Nothing,Float64}- Normalized axis fluxpsi_boundary_offset::Union{Nothing,Float64}- Boundary flux offsetpsi_axis_offset::Union{Nothing,Float64}- Axis flux offsetpsi_boundary_sign::Union{Nothing,Int}- Sign of boundary fluxpsi_axis_sign::Union{Nothing,Int}- Sign of axis fluxpsi_boundary_zero::Union{Nothing,Bool}- Whether boundary flux is zerormean::Union{Nothing,Float64}- Mean major radius [m]amean::Union{Nothing,Float64}- Mean minor radius [m]aratio::Union{Nothing,Float64}- Aspect ratio (R0/a)kappa::Union{Nothing,Float64}- Plasma elongationdelta1::Union{Nothing,Float64}- Upper triangularitydelta2::Union{Nothing,Float64}- Lower triangularitybt0::Union{Nothing,Float64}- Toroidal field at axis [T]crnt::Union{Nothing,Float64}- Plasma current [A]bwall::Union{Nothing,Float64}- Toroidal field at wall [T]verbose::Bool- Enable verbose outputdiagnose_src::Bool- Enable source data diagnosticsdiagnose_maxima::Bool- Enable extrema diagnosticsvolume::Union{Nothing,Float64}- Plasma volume [m³]betat::Union{Nothing,Float64}- Toroidal betabetan::Union{Nothing,Float64}- Normalized betabetaj::Union{Nothing,Float64}- Total betabetap1::Union{Nothing,Float64}- Poloidal beta (definition 1)betap2::Union{Nothing,Float64}- Poloidal beta (definition 2)betap3::Union{Nothing,Float64}- Poloidal beta (definition 3)li1::Union{Nothing,Float64}- Internal inductance (definition 1)li2::Union{Nothing,Float64}- Internal inductance (definition 2)li3::Union{Nothing,Float64}- Internal inductance (definition 3)
GeneralizedPerturbedEquilibrium.Equilibrium.FieldLineDerivParams — Type
FieldLineDerivParamsA struct to hold constant parameters for the ODE integration, making them easily accessible within the derivative function direct_fieldline_der!.
GeneralizedPerturbedEquilibrium.Equilibrium.GeometryProfileSplines — Type
GeometryProfileSplinesNamed 1D cubic spline interpolants for flux-surface-averaged geometric quantities. Built once during equilibrium construction so they are available to every downstream module without per-caller recomputation.
Fields
xs::Vector{Float64}: Shared ψ axis (matchesrzphi_xs)area_spline: Flux-surface area ∫ dA(ψ) [m²]avg_r_spline: Surface-average minor radius ⟨r⟩(ψ) [m]avg_R_spline: Surface-average major radius ⟨R⟩(ψ) [m]
GeneralizedPerturbedEquilibrium.Equilibrium.GeometryProfileSplines — Method
GeometryProfileSplines(xs, area_vals, avg_r_vals, avg_R_vals; extrap=ExtendExtrap())Construct GeometryProfileSplines from arrays of surface-averaged values defined on the shared ψ grid xs. Uses CubicFit boundary conditions to match ProfileSplines.
GeneralizedPerturbedEquilibrium.Equilibrium.InverseIngest — Type
InverseIngestThe serializable raw arrays and scalars captured by an inverse-equilibrium reader (read_chease_ascii, read_chease_binary) — everything needed to rebuild an InverseRunInput's splines without re-reading the original CHEASE file. Stored on InverseRunInput.ingest / PlasmaEquilibrium.ingest and reconstructed by build_inverse_from_ingest. See DirectIngest for the role this plays in the gpec.h5 rerun snapshot.
Fields
sq_xs::Vector{Float64}— normalized-ψ knots for the 1D profile splinesq_fs::Matrix{Float64}— 1D profile node valuesrz_xs::Vector{Float64}— ψ grid for the R, Z mapsrz_ys::Vector{Float64}— θ grid for the R, Z mapsR_nodes::Matrix{Float64}— R node values on the (ψ, θ) grid [m]Z_nodes::Matrix{Float64}— Z node values on the (ψ, θ) grid [m]ro::Float64— R of magnetic axis [m]zo::Float64— Z of magnetic axis [m]psio::Float64— total flux difference |ψaxis - ψboundary| [Wb/rad]
GeneralizedPerturbedEquilibrium.Equilibrium.InverseRunInput — Type
InverseRunInput(...)A container struct for inputs to the inverse_run function.
Fields
config::EquilibriumConfig- The equilibrium configuration objectsq_in::CubicSeriesInterpolant- 1D profile spline for F, P, qrz_in_xs::Vector{Float64}- ψ coordinate grid for rz_inrz_in_ys::Vector{Float64}- θ coordinate grid for rz_inrz_in_R::CubicInterpolantND- R coordinate interpolant [m]rz_in_Z::CubicInterpolantND- Z coordinate interpolant [m]ro::Float64- R-coordinate of magnetic axis [m]zo::Float64- Z-coordinate of magnetic axis [m]psio::Float64- Total flux difference |ψaxis - ψboundary| [Wb/rad]ingest::EquilibriumIngest- captured raw arrays for thegpec.h5rerun snapshot (anInverseIngestfor file-based reads, ornothingfor analytic equilibria)
GeneralizedPerturbedEquilibrium.Equilibrium.KineticProfileSplines — Type
KineticProfileSplinesNamed 1D cubic spline interpolants for kinetic profiles (densities, temperatures, ExB rotation, collisional diagnostics) loaded from an external kinetic.dat file. Mirrors the ProfileSplines pattern: each quantity is its own named spline plus a paired derivative view for the species the NTV kernel needs (densities and temperatures, used by wdian/wdiat in KineticForces/Torque.jl).
Fields
xs::Vector{Float64}: Shared ψ axis (the regular kinetic grid)ni_spline,ne_spline: Ion / electron number densities [m⁻³]Ti_spline,Te_spline: Ion / electron temperatures [J]omegaE_spline: ExB rotation ω_E [rad/s]loglam_spline: Coulomb logarithmnui_spline,nue_spline: Krook collision frequencies [s⁻¹]zeff_spline: Effective charge Z_effni_deriv,ne_deriv,Ti_deriv,Te_deriv: Derivative views for ψ-derivative access
GeneralizedPerturbedEquilibrium.Equilibrium.KineticProfileSplines — Method
KineticProfileSplines(xs, ni, ne, Ti, Te, omegaE, loglam, nui, nue, zeff;
extrap=ExtendExtrap())Build the kinetic profile splines from arrays of values defined on the shared ψ grid xs. Uses CubicFit boundary conditions to match ProfileSplines. Temperatures must already be in Joules.
GeneralizedPerturbedEquilibrium.Equilibrium.LargeAspectRatioConfig — Type
LargeAspectRatioConfig(...)A mutable struct holding parameters for the Large Aspect Ratio (LAR) plasma equilibrium model.
Fields:
lar_r0: The major radius of the plasma [m].lar_a: The minor radius of the plasma [m].beta0: The beta value on axis (normalized pressure).q0: The safety factor on axis.p_pres: The exponent for the pressure profile, defined asp00 * (1 - (r / a)^2)^p_pres.p_sig: The exponent that determines the shape of the current-related function profile.sigma_type: The type of sigma profile, can be "default" or "wesson". If "wesson", the sigma profile is defined assigma0 * (1 - (r / a)^2)^p_sig.mtau: The number of grid points in the poloidal direction.ma: The number of grid points in the radial direction.zeroth: If set to true, it neglects the Shafranov shift
GeneralizedPerturbedEquilibrium.Equilibrium.LargeAspectRatioConfig — Method
Build a LargeAspectRatioConfig from a parsed [LAR_INPUT] TOML table.
GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium — Type
PlasmaEquilibrium(...)The final, self-contained result of the equilibrium reconstruction. This object provides a complete representation of the processed plasma equilibrium in flux coordinates.
Fields
config::EquilibriumConfig: The equilibrium configuration object used for the reconstruction.params::EquilibriumParameters: Computed equilibrium parameters and diagnostics.profiles::ProfileSplines: Named 1D profile splines (F, P, dV/dψ, q) on normalized psi grid. Access values at grid points viaprofiles.F_spline.y[i], etc. Access derivatives viaprofiles.F_deriv.y[i]orprofiles.F_deriv(psi).geometry::GeometryProfileSplines: Named 1D splines for flux-surface-averaged geometry (area, ⟨r⟩, ⟨R⟩), populated automatically bycompute_geometry_profilesduring construction.Grid coordinates (shared by all rzphi/eqfun interpolants):
rzphi_xs::Vector{Float64}: ψ coordinates (length mpsi+1)rzphi_ys::Vector{Float64}: θ coordinates (length mtheta+1)
Geometric quantities (rzphi, 4 interpolants): 2D cubic interpolants for flux-coordinate mapping with periodic BC in theta.
- x value: normalized ψ
- y value: SFL poloidal angle ∈ [0, 1]
rzphi_rsquared::CubicInterpolantND: r_coord² = (R - ro)² + (Z - zo)²rzphi_offset::CubicInterpolantND: η/(2π) - θₙₑw (angle offset)rzphi_nu::CubicInterpolantND: ν in ϕ = 2πζ + ν(ψ, θ)rzphi_jac::CubicInterpolantND: Jacobian
Physics quantities (eqfun, 3 interpolants): 2D cubic interpolants storing local physics and geometric quantities.
- x value: normalized ψ
- y value: SFL poloidal angle θₙₑw
eqfun_B::CubicInterpolantND: Total magnetic field strength [T]eqfun_metric1::CubicInterpolantND: (e₁⋅e₂ + q⋅e₃⋅e₁)/(J⋅B²)eqfun_metric2::CubicInterpolantND: (e₂⋅e₃ + q⋅e₃⋅e₃)/(J⋅B²)
ro::Float64: R-coordinate of the magnetic axis [m]zo::Float64: Z-coordinate of the magnetic axis [m]psio::Float64: Total flux difference |Ψaxis - Ψboundary| [Weber/radian]ingest::EquilibriumIngest: raw arrays forwarded from the equilibrium input for thegpec.h5rerun snapshot — aDirectIngest/InverseIngestfor file-based equilibria, ornothingfor analytic ones (regenerated from their TOML section on replay)
GeneralizedPerturbedEquilibrium.Equilibrium.ProfileSplines — Type
ProfileSplinesNamed 1D cubic spline interpolants for equilibrium profiles. Each profile is stored as a separate spline for code clarity.
Fields
xs::Vector{Float64}: Shared x-axis (normalized psi)F_spline: 2π*F (toroidal flux function, where F = R * B_toroidal)P_spline: μ₀*P (plasma pressure × μ₀)dVdpsi_spline: dV/dψ (volume derivative)q_spline: q (safety factor)
Derivative Interpolants (for continuous derivative evaluation)
F_deriv,P_deriv,dVdpsi_deriv,q_deriv: First derivative interpolants
Notes
- Node values at grid points: Access via
spline.y[i] - Derivative at any point: Call
deriv(x)(derivative views are callable) - Grid: Access via
xsfield orspline.cache.x
GeneralizedPerturbedEquilibrium.Equilibrium.ProfileSplines — Method
ProfileSplines(xs, F_vals, P_vals, dVdpsi_vals, q_vals; extrap=ExtendExtrap())Create ProfileSplines from arrays of profile values. Uses CubicFit boundary conditions with extension extrapolation.
GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig — Type
SolovevConfig(...)A mutable struct holding parameters for the Solev'ev (SOL) plasma equilibrium model.
Fields:
mr: number of radial grid zonesmz: number of axial grid zonesma: number of flux grid zonese: elongationa: minor radiusr0: major radiusq0: safety factor at the o-pointp0fac: scale on-axis pressure (P-> P+P0*p0fac. beta changes. Phi,q constant)b0fac: scale toroidal field at constant beta (sPhi,sf,s^2*P. bt changes. Shape,beta constant)f0fac: scale toroidal field at constant pressure (s*f. beta,q changes. Phi,p,bp constant)
GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig — Method
Build a SolovevConfig from a parsed [SOL_INPUT] TOML table.
GeneralizedPerturbedEquilibrium.Equilibrium.TJAnalyticConfig — Type
TJAnalyticConfig(...)Parameters for the TJ-analytic cylindrical large-aspect-ratio equilibrium model — a GPEC adaptation of the analytic profile family used by R. Fitzpatrick's TJ code (https://github.com/rfitzp/TJ). We follow the same analytic-profile parameterization (ψ-ODE in dimensionless r/a, f₁ for q, power-law pressure) for the inner cylindrical core and connect it to GPEC's direct-GS pipeline; this is NOT a re-implementation of TJ.
The model uses analytic profiles with exact control of both the on-axis and edge safety factors. The q profile is determined by:
f1(r) = [1 - (1-r²)^ν] / (ν·qc)
q(r) = r² / f1(r)where ν = qa/qc is the current peaking parameter, qc is the axis q, and qa is the edge q. All lengths are normalized to R₀, fields to B₀. The pressure profile is p₂(r) = pc·(1-r²)^μ.
Reference: R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ
GeneralizedPerturbedEquilibrium.Equilibrium.TJAnalyticConfig — Method
Build a TJAnalyticConfig from a parsed [TJ_ANALYTIC_INPUT] TOML table.
GeneralizedPerturbedEquilibrium.Equilibrium.TJAnalyticShapeParams — Type
Internal parameter bundle for the TJ-analytic shape ODE (ψ, g₂, H₁, H₁', f₃) — GPEC adaptation of the analytic shape ODE used in R. Fitzpatrick's TJ code (https://github.com/rfitzp/TJ). Built once per tj_analytic_run / tj_analytic_run_direct call so both pipelines share identical numerics.
Fields:
- physical: a, R0, qc, mu, pc, B0
- derived: epsa2 = (a/R0)²
- near-axis BC constants: rmin, x0 = rmin, r0 = rmin·a, f1c = 1/qc, p2ppc = d²p₂/dx²|_0 = −2·μ·pc
GeneralizedPerturbedEquilibrium.Equilibrium._build_psi_grid — Method
_build_psi_grid(equil_params, psilow, psihigh, fieldline_int, raw_profile, ro, zo, rs2)Resolve mpsi and build psi_nodes for any supported grid_type.
For "log_asymptotic" with mpsi=0, estimates the separatrix log slope from two probe integrations and computes the minimum knot count for psi_accuracy. For "ldp", uses the sin²-spaced grid. Shared by direct_fieldline_int and efit_by_inversion solvers.
GeneralizedPerturbedEquilibrium.Equilibrium._cubic_resample — Method
Cubic-spline resample matching Fortran pentrc/inputs.f90:215-232: build a cubic spline on the (irregular) input grid, then evaluate at the regular psinew grid. Out-of-range points use ExtendExtrap (smooth cubic extrapolation), matching Fortran's `splinefit(...,"extrap")`.
GeneralizedPerturbedEquilibrium.Equilibrium._estimate_log_slope — Method
_estimate_log_slope(fieldline_int, raw_profile, ro, zo, rs2, psihigh)Estimate the logarithmic slope A from two field-line probe integrations near the separatrix, using the asymptotic form q(ψ) ≃ −A·ln(1−ψ) (Fitzpatrick 2024, eq. 19).
Returns A = |Δq| / ln(2). Falls back to A=2.0 on integration failure.
GeneralizedPerturbedEquilibrium.Equilibrium._estimate_mid_spacing — Method
_estimate_mid_spacing(sq_in, psi_split_core, psi_split_edge, tau)Estimate the uniform knot spacing needed in the middle ψ region to resolve the sharpest profile feature (usually the pressure pedestal) to relative accuracy tau.
Uses central-difference second derivatives of all 4 sqin profiles on a 300-point sample. The required spacing for profile k is hk = sqrt(8τ / d2normk) where d2normk is max|f''| / max|f| over [psisplitcore, psisplitedge].
GeneralizedPerturbedEquilibrium.Equilibrium._is_closed_curve — Method
_is_closed_curve(curve)Returns true if the Contour.jl curve is closed (first ≈ last vertex).
GeneralizedPerturbedEquilibrium.Equilibrium._read_1d_gfile_format — Method
read1dgfileformat(linesblock, numvalues)
Internal helper function to parse Fortran-style fixed-width numerical blocks from a vector of strings.
Arguments:
lines_block: AVector{String}containing the lines to parse.num_values: The total number ofFloat64values to read from the block.
Returns:
- A
Vector{Float64}containing the parsed values.
GeneralizedPerturbedEquilibrium.Equilibrium._read_kinetic_table — Method
Internal helper: parse the raw 6-column kinetic profile table from disk, filtering out non-numeric header rows. Returns six independent column views. # comment lines (e.g. a provenance header) are stripped so they cannot widen the parsed matrix and pad the data rows.
GeneralizedPerturbedEquilibrium.Equilibrium._winding_number_verts — Method
_winding_number_verts(verts, r_test, z_test)Winding number point-in-polygon test (Sunday algorithm) on a Ctr.vertices vector. Returns the signed winding count; nonzero means the test point is inside the polygon. More robust than ray-casting for near-collinear edges (e.g., surfaces near an x-point).
GeneralizedPerturbedEquilibrium.Equilibrium.adaptive_grid_params — Method
adaptive_grid_params(raw_profile, ro, zo, r_lo, r_hi, z_lo, z_hi,
psilow, Δψ, bbox_curve) → (nr, nz, β_r, β_z)Physics-based Cartesian grid sizing for contour-tracing equilibrium reconstruction.
Computes required cell widths from the bilinear interpolation error bound δψ ≈ (dR²·|ψRR| + dZ²·|ψZZ|)/8 ≤ Δψ → dRreq = √(8Δψ/|ψRR|) sampled at each LCFS vertex (reusing bbox_curve), plus an axis constraint ensuring the innermost surface (psilow) has ≥ 5 cells per semi-axis on the global grid.
β = acosh(hmax/hmin) from the ratio of coarsest to finest required cell widths. n = 1 + ⌈span·sinch(β)/h_min⌉ where sinch(β) = β/sinh(β).
GeneralizedPerturbedEquilibrium.Equilibrium.arclength_fieldline_der! — Method
arclength_fieldline_der!(dy, y, params, s)RHS for the arc-length-parameterized level-set ODE.
State y = [R, Z, ∫dl/Bp, ∫dl/(R²Bp), ∫jac·dl/Bp]. The independent variable s is arc length [m].
The tangent direction (dR/ds, dZ/ds) = (∂ψ/∂Z, −∂ψ/∂R) / |∇ψ| follows the ψ = const level set counterclockwise, staying on the flux surface without correction.
Near x-points where Bp → 0, the position integration (dy[1:2]) remains well-behaved because grad_norm never appears in the denominator alone. dy[3:5] use abstol=1e20 so the solver never restricts step size for them.
GeneralizedPerturbedEquilibrium.Equilibrium.arclength_fieldline_int — Method
arclength_fieldline_int(psifac, raw_profile, ro, zo, rs2)Arc-length-parameterized flux surface integration. Drop-in replacement for direct_fieldline_int with identical return format:
y_out[:, 1]: geometric angle η ∈ 0 to 2π (CCW from outboard midplane)y_out[:, 2]: accumulated ∫dl/Bpy_out[:, 3]: rfac = √((R−ro)² + (Z−zo)²)y_out[:, 4]: accumulated ∫dl/(R²Bp)y_out[:, 5]: accumulated ∫jac·dl/Bp
The ODE is terminated by a ContinuousCallback that detects the return to the outboard midplane (Z = zo, R > ro) after a minimum arc-length guard.
GeneralizedPerturbedEquilibrium.Equilibrium.area_to_rootarea_weight — Method
area_to_rootarea_weight(equil, psi, ft) -> Matrix{ComplexF64}Inverse of rootarea_to_area_weight: the area-weighted → root-area-weighted field operator √A·Σ⁻¹ = √jarea · inv(sqrtamat), mapping b̄ → b̃ (b̃ = √A·Σ⁻¹·b̄). [Pharr 2026]
GeneralizedPerturbedEquilibrium.Equilibrium.build_direct_from_ingest — Method
build_direct_from_ingest(config::EquilibriumConfig, ingest::DirectIngest) -> DirectRunInputRebuild a DirectRunInput from a DirectIngest captured by read_efit/read_imas (or restored from input/raw_inputs/equilibrium/ inside gpec.h5). Inverse of that capture: reconstructs the splines so the rerun path skips the g-file/IMAS parse, reusing the existing solver dispatch.
GeneralizedPerturbedEquilibrium.Equilibrium.build_inverse_from_ingest — Method
build_inverse_from_ingest(config::EquilibriumConfig, ingest::InverseIngest) -> InverseRunInputRebuild an InverseRunInput from an InverseIngest captured by read_chease_ascii/read_chease_binary. Inverse of that capture; used by the rerun path to skip the CHEASE parse.
GeneralizedPerturbedEquilibrium.Equilibrium.clamp_psihigh_to_separatrix — Method
clamp_psihigh_to_separatrix(raw_profile) -> (clamped_psihigh, was_adjusted)Binary-searches for the highest psihigh ≤ raw_profile.config.psihigh at which the ψ level set is still a closed curve in the EFIT grid. Returns the safe value and a Bool indicating whether any clamping occurred.
GeneralizedPerturbedEquilibrium.Equilibrium.classify_topology — Method
classify_topology(raw_profile, psio; xpt_threshold=0.05) → SymbolClassify the plasma topology as :limited, :sn_lower, :sn_upper, or :double_null by evaluating ψ on the (R, Z) grid and finding the minimum in each Z half-domain.
An x-point is detected when the minimum ψ in a half-domain falls below xpt_threshold * psio.
GeneralizedPerturbedEquilibrium.Equilibrium.compute_geometry_profiles — Method
compute_geometry_profiles(pe::PlasmaEquilibrium) → GeometryProfileSplinesConvenience method that pulls the rzphi grids and interpolants directly off a PlasmaEquilibrium. Useful for tests and after-the-fact recomputation.
GeneralizedPerturbedEquilibrium.Equilibrium.compute_geometry_profiles — Method
compute_geometry_profiles(rzphi_xs, rzphi_ys,
rzphi_rsquared, rzphi_offset, rzphi_jac, ro)
→ GeometryProfileSplinesBuild named cubic splines for flux-surface-averaged area, ⟨r⟩, and ⟨R⟩ from the equilibrium's 2D rzphi interpolants. Uses trapezoidal sums in θ across the non-periodic θ nodes (the trailing duplicate at θ = 1 is skipped, mthsurf = length(rzphi_ys) - 1).
The integrand delpsi = |∇ψ| (in geometry units) is built from the same nodal derivatives the GS-residual diagnostic uses; see Equilibrium.jl:407-419 and the Fortran reference at pentrc/dcon_interface.f:1199-1254.
GeneralizedPerturbedEquilibrium.Equilibrium.compute_sqrt_jac_delpsi — Method
compute_sqrt_jac_delpsi(equil, psi, mtheta) -> Vector{Float64}Compute √(J·|∇ψ|) at mtheta equally-spaced θ points on the flux surface at psi. This is the √weight function that maps a field component b(θ) to its root-area-weighted form √(J|∇ψ|)·b(θ) in θ-space.
GeneralizedPerturbedEquilibrium.Equilibrium.compute_sqrtamat — Method
compute_sqrtamat(equil, psi, ft) -> Matrix{ComplexF64}Build the √A convolution matrix sqrtamat. Produces a Hermitian Toeplitz matrix with entries sqrtamat[m',k] = ŵ{mk − m'} where ŵn = (1/N)Σ wj exp(+inθ_j) and w(θ) = √(J·|∇ψ|).
Operationally, sqrtamat is the mode-space √weight operator: for a field b with Fourier coefficients bfft, it satisfies the identity `‖sqrtamat·bfft‖² = N² · ∫ |b|² · J|∇ψ| dθwhich is Jacobian-invariant on a given flux surface (seescripts/testpowernorm_invariance.jl`).
The backward step uses exp(−imθ)/(1/N) normalization paired with the Julia forward FT exp(+imθ) so that round-trip = identity and the convolution structure is correct.
GeneralizedPerturbedEquilibrium.Equilibrium.direct_fieldline_der! — Method
direct_fieldline_der!(dy, y, params, eta)The derivative function for the field-line integration ODE. This is passed to the DifferentialEquations.jl solver. This is a Julia adaptation of the Fortran direct_fl_der subroutine, with an added safeguard against division by zero.
Arguments:
dy: The derivative vector (output, modified in-place).y: The state vector[∫(dl/Bp), rfac, ∫(dl/(R²Bp)), ∫(jac*dl/Bp)].params: AFieldLineDerivParamsstruct with all necessary parameters.eta: The independent variable (geometric angleη).
GeneralizedPerturbedEquilibrium.Equilibrium.direct_fieldline_int — Method
direct_fieldline_int(psifac, raw_profile, ro, zo, rs2)Performs the field-line integration for a single flux surface. This is a Julia adaptation of the Fortran direct_fl_int subroutine. Note that the array y_out is now indexed from 1:5 rather than 0:4 as in Fortran.
Arguments:
psifac: normalized psi value for the surface (ψ_norm).raw_profile:DirectRunInputobject containing splines and parameters.ro,zo: Coordinates of the magnetic axis [m].rs2: R-coordinate of the outboard separatrix [m].
Returns:
- `y_out`: A matrix containing the integrated quantities vs. the geometric angle `η`.
- `y_out[:, 1]`: η (geometric poloidal angle)
- `y_out[:, 2]`: ∫(dl/Bp)
- `y_out[:, 3]`: rfac (radial distance from magnetic axis)
- `y_out[:, 4]`: ∫(dl/(R²Bp))
- `y_out[:, 5]`: ∫(jac*dl/Bp)bfield: ADirectBFieldobject with values at the integration start point.
GeneralizedPerturbedEquilibrium.Equilibrium.direct_get_bfield! — Method
direct_get_bfield!(bf_out, r, z, psi_in, sq_in, sq_in_deriv, psio; derivs=0)Calculates the magnetic field and its derivatives at a given (R,Z) point. The results are stored in-place in the bf_out object. This is equivalent to the direct_get_bfield subroutine in the Fortran code, adapted for the Julia spline implementation.
Arguments:
bf_out: A mutableDirectBFieldstruct to store the resultsr: R-coordinate to evaluate atz: Z-coordinate to evaluate atpsi_in: 2D cubic interpolant for poloidal fluxψ(R,Z)sq_in: 1D cubic spline for profilesF(ψ_norm)andP(ψ_norm)sq_in_deriv: Pre-computed derivative view of sq_inpsio: total toroidal fluxderivs: An integer specifying number of derivatives to compute (0, 1, or 2)
GeneralizedPerturbedEquilibrium.Equilibrium.direct_position! — Method
direct_position!(raw_profile)Finds the key geometric locations of the equilibrium: the magnetic axis (O-point) and the inboard/outboard separatrix crossings on the midplane. It also updates the spline representing the poloidal flux ψ(R,Z) based on the new magnetic axis location. This function performs the same overall function as the Fortran direct_position subroutine with better iteration control and error handling. We have also added a helper function for separatrix finding.
Arguments:
raw_profile: ADirectRunInputobject containing splines and parameters.
Returns:
ro: R-coordinate of the magnetic axis [m].zo: Z-coordinate of the magnetic axis [m].rs1: R-coordinate of the inboard separatrix crossing [m].rs2: R-coordinate of the outboard separatrix crossing [m].psi_in_new: returns psi_in renormalized by * psio/psi(ro,zo)
GeneralizedPerturbedEquilibrium.Equilibrium.direct_refine — Method
direct_refine(rfac, eta, psi0, params)Refines the radial distance rfac at a given angle eta to ensure the point lies exactly on the target flux surface psi0.
Arguments:
rfac: The current guess for the radial distance from the magnetic axis.eta: The geometric poloidal angle.psi0: The targetψvalue for the flux surface.params: AFieldLineDerivParamsstruct.
Returns:
- The refined
rfacvalue.
GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_global_parameters! — Method
equilibrium_global_parameters!(pe::PlasmaEquilibrium)Computes and populates global equilibrium parameters in the PlasmaEquilibrium struct, such as rmean, amean, kappa, bt0, crnt, betat, betan, li1, etc. Performs the same function as equiloutglobal in the Fortran code.
GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_gse! — Method
equilibrium_gse!(equil::PlasmaEquilibrium)Diagnoses the Grad-Shafranov solution by computing the residual of the Grad-Shafranov equation across the grid and writing diagnostic data to HDF5 files. Performs the same function as equiloutgse in the Fortran code.
GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_qfind! — Method
equilibrium_qfind!(equil::PlasmaEquilibrium)Finds the extrema of the safety factor profile q(ψ) in the plasma equilibrium and computes derived q-values such as q0, qmin, qmax, qa, and q95. Performs the same function as equiloutqfind in the Fortran code.
GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_separatrix_find! — Method
equilibrium_separatrix_find!(pe::PlasmaEquilibrium)Finds the separatrix locations in the plasma equilibrium (rsep, zsep, rext, zext). Performs the same function as equiloutsep_find in the Fortran code.
GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_solver — Function
equilibrium_solver(raw_profile)The main driver for the direct equilibrium reconstruction. It orchestrates the entire process from finding the magnetic axis to integrating along field lines and constructing the final coordinate and physics quantity splines. This performs the same overall function as the Fortran direct_run subroutine, with better checks for numerical robustness.
Arguments:
raw_profile: ADirectRunInputobject containing the initial splines (psi_in,sq_in) and run parameters (equil_input).
Returns:
- A
PlasmaEquilibriumobject containing the final, processed equilibrium data, including the profile spline (sq), the coordinate mapping spline (rzphi), and the physics quantity spline (eqfun).
GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_solver_by_inversion — Method
equilibrium_solver_by_inversion(raw_profile; resolution_factor=1.0,
refine=nothing, β_r=nothing, β_z=nothing)Driver for contour-tracing equilibrium reconstruction.
Grid parameters (n and β in each direction) are computed adaptively from the bilinear interpolation error bound along the LCFS and an axis constraint. resolution_factor scales both n values uniformly (higher → more points, same β).
Legacy keyword arguments refine, β_r, β_z override the adaptive values when provided (for benchmark sweeps and backward compatibility).
Select via eq_type = "efit_by_inversion" in gpec.toml.
GeneralizedPerturbedEquilibrium.Equilibrium.flux_surface_area — Method
flux_surface_area(equil, psi, mtheta) -> Float64Compute ∫ J·|∇ψ| dθ over the flux surface at psi using mtheta equally-spaced θ points and the trapezoidal rule with periodic end correction.
GeneralizedPerturbedEquilibrium.Equilibrium.flux_surface_metric — Method
flux_surface_metric(equil, psi, theta; hint=(Ref(1), Ref(1))) -> (; r, jac, delpsi)Evaluate the flux-surface geometry at (psi, theta): major radius r, Jacobian jac, and the flux-gradient magnitude delpsi = |∇ψ|. The metric components w11/w12 follow GPEC's gpout.f convention for the rzphi straight-field-line coordinates. This is the single source of truth for the |∇ψ| computation reused across the FFS and PE modules.
GeneralizedPerturbedEquilibrium.Equilibrium.inverse_extrap — Method
inverse_extrap(xx::Matrix{Float64}, ff::Matrix{Float64}, x::Float64) -> Vector{Float64}Performs component-wise Lagrange extrapolation for a vector-valued function.
Arguments:
xx: A (m × n) matrix where each row contains the x-values for each component.ff: A (m × n) matrix where each row contains function values at the correspondingxx.x: A scalar Float64 value at which to extrapolate.
Returns:
- A vector of length n representing the extrapolated function values at
x.
GeneralizedPerturbedEquilibrium.Equilibrium.lar_init_conditions — Method
lar_init_conditions(rmin, sigma_type, params)Initializes the starting radius and state vector for solving the LAR ODE system. Also evaluates the initial derivative using the analytic model.
Arguments:
rmin: Normalized starting radius (as a fraction oflar_a).lar_input: ALargeAspectRatioConfigobject containing equilibrium parameters.
Returns:
r: Physical radius corresponding tormin * lar_a.y: Initial state vector of length 5.
GeneralizedPerturbedEquilibrium.Equilibrium.load_kinetic_profiles — Method
load_kinetic_profiles(kinetic_file::AbstractString;
zi::Int=1, zimp::Int=6, mi::Int=2, mimp::Int=12,
density_factor::Float64=1.0, temperature_factor::Float64=1.0,
ExB_rotation_factor::Float64=1.0, toroidal_rotation_factor::Float64=1.0,
chi1::Union{Nothing,Float64}=nothing)
→ KineticProfileSplinesParse an ASCII kinetic profile file, interpolate onto a regular 101-point ψ grid, optionally apply profile scaling knobs, derive collisional / Z_eff diagnostics, and return a KineticProfileSplines with independent named cubic splines.
Expected file format
Six whitespace-separated columns (header rows are filtered out):
psi_n n_i[m^-3] n_e[m^-3] T_i[eV] T_e[eV] omega_E[rad/s]Arguments
kinetic_file: Path to the ASCII kinetic profile filezi,zimp: Main ion and impurity charge numbersmi,mimp: Main ion and impurity mass numbers (in proton masses)density_factor: Density scaling factor (applied to ni, ne)temperature_factor: Temperature scaling factor (applied to Ti, Te)ExB_rotation_factor: ExB rotation scaling factor (applied to omegaE after rotation reform)toroidal_rotation_factor: Toroidal rotation scaling factor (scales total wphi = omegaE + wdian + wdiat)chi1: Poloidal flux normalization2π·ψ₀— required whendensity_factor,temperature_factor, ortoroidal_rotation_factordiffer from 1.0
Scaling sequence
When any of density_factor, temperature_factor, toroidal_rotation_factor differ from 1.0:
- Build first-pass cubic splines from unscaled profiles (for derivatives)
- Compute diamagnetic frequencies
wdian,wdiatand total toroidal rotationwphi - Scale:
wdian_new = temperature_factor * wdian,wdiat_new = temperature_factor * wdiat(density_factorcancels inT*(dn/dψ)/n) - Reform:
omegaE = toroidal_rotation_factor * wphi - wdian_new - wdiat_new - Scale density/temperature arrays:
ni *= density_factor,Ti *= temperature_factor, etc.
Then ExB_rotation_factor is applied independently: omegaE *= ExB_rotation_factor.
Collisionality is recomputed from the (possibly scaled) profiles. This differs from Fortran PENTRC, which computes collisionality from unscaled profiles. Use nufac (in KineticForcesControl) for independent collisionality scaling.
GeneralizedPerturbedEquilibrium.Equilibrium.make_multi_stretched_grid — Method
make_multi_stretched_grid(lo, hi, centers, n, β) → Vector{Float64}Build a 1D grid of length n on [lo, hi] with sinh-stretching toward each point in centers. The domain is split at each interior concentration point; each resulting sub-interval uses:
- symmetric (double-ended) sinh when both endpoints are concentration points
- one-sided sinh fine at the concentration-point end otherwise
Domain boundaries (lo, hi) are treated as concentration points only if they appear in centers. This produces a grid that is simultaneously fine near the magnetic axis (an interior concentration point) and near x-points (which lie at domain boundaries for the clipped plasma bounding box). Returns a uniform grid if β ≤ 0.
GeneralizedPerturbedEquilibrium.Equilibrium.make_optimal_mpsi — Method
make_optimal_mpsi(psilow, psihigh, A, sq_in; tau, psi_split_core, psi_split_edge)Compute the minimum number of radial knots needed to achieve target accuracy τ in q, given the separatrix log slope A. Three-region geometric grid: core, pedestal, far edge. The middle-region spacing is driven by profile curvature (P, F, dV/dψ, q) via sq_in.
GeneralizedPerturbedEquilibrium.Equilibrium.make_optimal_psi_grid — Method
make_optimal_psi_grid(psilow, psihigh, N_core, N_mid, N_edge; psi_split_core, psi_split_edge)Build a three-region ψ grid with (Ncore + Nmid + N_edge)+1 knots, using the counts directly as provided — core and edge are geometric in log(ψ) and log(1−ψ) respectively; middle is uniform in ψ.
GeneralizedPerturbedEquilibrium.Equilibrium.make_stretched_r_grid — Method
make_stretched_r_grid(r_lo, r_hi, ro, nr, β_r) → Vector{Float64}Grid of nr points on [r_lo, r_hi] concentrated near ro (magnetic axis). Thin wrapper around make_multi_stretched_grid.
GeneralizedPerturbedEquilibrium.Equilibrium.make_stretched_z_grid — Method
make_stretched_z_grid(z_lo, z_hi, zo, nz, topology, β_z) → Vector{Float64}Grid of nz points on [z_lo, z_hi] concentrated near zo (magnetic axis) and near each x-point boundary (at z_lo for :sn_lower/:double_null, at z_hi for :sn_upper/:double_null). Thin wrapper around make_multi_stretched_grid.
GeneralizedPerturbedEquilibrium.Equilibrium.read_chease_ascii — Method
read_chease_ascii(config)Parses a ascii CHEASE file, creates initial 1D and 2D splines, finds magnetic axis, and bundles them into a InverseRunInput object.
Arguments:
config: TheEquilibriumConfigobject containing the filename and parameters.
Returns:
- A
InverseRunInputobject ready for the inverse solver.
GeneralizedPerturbedEquilibrium.Equilibrium.read_chease_binary — Method
read_chease_binary(equil_config)Parses a binary CHEASE file, creates initial 1D and 2D splines with proper normalization (R0, B0 scaling), and bundles them into a InverseRunInput object.
GeneralizedPerturbedEquilibrium.Equilibrium.read_efit — Method
_read_efit(equil_in)Parses an EFIT g-file, creates initial 1D and 2D splines, and bundles them into a DirectRunInput object.
Arguments:
equil_in: TheEquilInputobject containing the filename and parameters.
Returns:
- A
DirectRunInputobject ready for the direct solver.
GeneralizedPerturbedEquilibrium.Equilibrium.read_imas — Method
read_imas(config::EquilibriumConfig, dd)Load an equilibrium from an IMAS data dictionary and return a DirectRunInput.
The dd.equilibrium.time_slice[] is used (active time slice). Poloidal flux is converted from the IMAS COCOS convention (set by config.imas_cocos) to the internal COCOS 2 convention:
imas_cocos = 11(default, IMAS standard): divide ψ by 2πimas_cocos = 2(GPEC internal): no conversion
Arguments
config:EquilibriumConfigwitheq_type = "imas"andimas_cocosset.dd: populatedIMASdd.ddwithdd.equilibrium.time_slice[]containing:global_quantities.psi_axis,global_quantities.psi_boundaryprofiles_1d.psi,profiles_1d.f,profiles_1d.pressure,profiles_1d.qprofiles_2d[1].grid.dim1(R),profiles_2d[1].grid.dim2(Z),profiles_2d[1].psi
GeneralizedPerturbedEquilibrium.Equilibrium.resample_contour_to_theta_grid! — Method
resample_contour_to_theta_grid!(R_out, Z_out, curve, ro, zo, theta_grid)Resample a Contour.jl curve onto a uniform geometric-angle grid, writing results directly into the pre-allocated R_out and Z_out vectors (views into Rtable/Ztable).
Uses a polar representation: fits a cubic spline on log ρ(η) where ρi = sqrt((Ri − ro)² + (Zi − zo)²) and ηi = atan2(Zi − zo, Ri − ro). This guarantees exp(logρ_spl) > 0 everywhere — a self-intersecting contour is impossible by construction — making it robust near x-points where R(η) and Z(η) splines tend to overshoot between widely-spaced knots.
GeneralizedPerturbedEquilibrium.Equilibrium.rootarea_to_area_weight — Method
rootarea_to_area_weight(equil, psi, ft) -> Matrix{ComplexF64}Build the root-area-weighted → area-weighted field operator Σ/√A = sqrtamat ./ √jarea at the flux surface psi, where jarea = ∫ J|∇ψ| dθ is the scalar flux-surface area. It maps the coordinate-invariant root-area-weighted field b̃ to the area-weighted field b̄: b̄ = (Σ/√A)·b̃ (both in tesla). Poloidal flux is the scalar product Φ = A·b̄ (so Φ = Σ·√A·b̃), recovered only when a user supplies/requests flux — it is never stored.
The √-weight (b̃) basis is the one in which operator singular values / spectra are independent of the straight-field-line (working) coordinate — see field_space_response_matrices. b̄ is not coordinate-invariant; it is only a field/flux recovery view. [Pharr 2026]
GeneralizedPerturbedEquilibrium.Equilibrium.select_plasma_contour — Method
select_plasma_contour(curve_list, ro, zo)From the list of Contour.jl curves at a given ψ level, return the closed curve whose interior contains the magnetic axis (ro, zo).
Returns nothing if no qualifying closed curve is found.
GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium — Function
setup_equilibrium(eq_config::EquilibriumConfig)Read an equilibrium file, run the appropriate solver, and return the processed PlasmaEquilibrium with global parameters, q-profile, and GSE diagnostics.
GeneralizedPerturbedEquilibrium.Equilibrium.sol_run — Method
This function handles the Solovev analytical equilibrium model, transforming the input parameters into the necessary splines and scalar values for equilibrium construction. This is a Julia version of the Fortran code in sol.f, with no major differences except for arrays going from 0:n to 1:n+1.
Arguments:
mr: Number of radial grid zonesmz: Number of axial grid zonesma: Number of flux grid zonese: Elongationa: Minor radiusr0: Major radiusq0: Safety factor at the o-pointp0fac: Scales on axis pressure (s*P. beta changes. Phi,q constant)b0fac: Scales on toroidal field (sPhi,sf,s^2*P. bt changes. Shape,beta constant)f0fac: Scales on toroidal field (s*f. bt,q changes. Phi,p,bp constant)
Returns:
DirectRunInputobject
GeneralizedPerturbedEquilibrium.Equilibrium.tj_analytic_f1 — Method
tj_analytic_f1(x, nu, qc)TJ-analytic poloidal flux function f1(x) where x = r/a, following the analytic-profile parameterization of R. Fitzpatrick's TJ code (https://github.com/rfitzp/TJ). Uses a Taylor expansion near the axis for numerical stability.
Reference: R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ
GeneralizedPerturbedEquilibrium.Equilibrium.tj_analytic_f1p — Method
tj_analytic_f1p(x, nu, qc)Derivative of the TJ-analytic f1 with respect to x (= r/a). See R. Fitzpatrick's TJ code (https://github.com/rfitzp/TJ) for the original parameterization.
GeneralizedPerturbedEquilibrium.Equilibrium.tj_analytic_find_nu — Method
TJ-analytic ν root-find (Fitzpatrick's Setnu / GetNu in https://github.com/rfitzp/TJ): solve for ν so that q₂(x=1) matches qa_target.
q₂ = x²·(1+εa²·g₂)·exp(−εa²·f3/f1)/f1; at x=1 and low β this picks up an O(εa²) correction relative to the lowest-order guess ν = qa/qc, which matters for the TJ-analytic benchmark at large ε. Falls back to the lowest-order ν if the bracket search diverges.
GeneralizedPerturbedEquilibrium.Equilibrium.tj_analytic_run — Method
tj_analytic_run(equil_input, tj_input)Construct a cylindrical tokamak equilibrium using the TJ-analytic model — GPEC's adaptation of the analytic-profile family used in R. Fitzpatrick's TJ code (https://github.com/rfitzp/TJ).
Profiles are analytic:
f1(x) = [1 - (1-x²)^ν] / (ν·qc), p2(x) = pc·(1-x²)^μ, x = r/awith ν = qa/qc. The 2D geometry is built from the TJ-analytic inverse aspect-ratio expansion. With zero edge shaping (Hna = Vna = 0) — the TJ-analytic benchmark configuration of Fitzpatrick's TJ — flux surfaces are shifted circles
R(r,θ) = R₀ + Δ(r) + α(r)·r·cos θ
Z(r,θ) = α(r)·r·sin θwhere Δ and α come from the shaping ODE for (g₂, H₁, H₁') (same equations as Fitzpatrick's TJ shape ODE):
Δ(r) = R₀·εa²·H₁(x) (Shafranov shift)
α(r) = 1 − εa²·(x²/8 − H₁/2) (from L(x) = x³/8 − x·H₁/2)
εa = a/R₀The higher-order toroidal-flux correction g₂ enters the output F profile as F = R₀·B₀·(1 + εa²·g₂), and the higher-order poloidal flux f₃ enters the safety factor as q₂ = x²·(1 + εa²·g₂)·exp(−εa²·f₃/f1)/f1.
The (n ≥ 2) horizontal/vertical shaping harmonics Hₙ(r), Vₙ(r) are not yet included; they are zero in the TJ-analytic benchmark scans.
Reference: R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ
GeneralizedPerturbedEquilibrium.Equilibrium.tj_analytic_run_direct — Method
tj_analytic_run_direct(equil_input, tj_input; nrbox=257, nzbox=257, rc=1.2)Option B pipeline: construct ψ(R, Z) on a 2D grid from the TJ-analytic model — GPEC's adaptation of R. Fitzpatrick's TJ code analytic-profile family (https://github.com/rfitzp/TJ) — and return a DirectRunInput so the equilibrium is processed by the direct-GS solver (same path as the geqdsk-based scans).
Using the inverse pipeline on just the first-order Shafranov-shifted-circle geometry systematically under-drives the external kink at large ε because the inverse solver consumes the prescribed q₂ profile and never recomputes q from geometry. The direct pipeline, in contrast, line-integrates F·∮dθ/(R²·Bp) on the 2D ψ(R,Z) field, so higher-order geometric effects (buried in the shape of ψ away from the axis) feed back into q and δW. Reproducing the full geqdsk-equivalent path therefore requires rebuilding ψ(R,Z) from the analytic model itself — not just the flux-surface coordinates — including the vacuum region outside the plasma.
The benchmark keeps edge shaping Hna = Vna = 0, so the ODE-integrated shape harmonics Hₙ, Vₙ for n ≥ 2 are rescaled to zero; only the H₁ Shafranov shift contributes. ψ(R, Z) is constructed by:
- for each grid point, iterating the map (R, Z) → (r, w) 10× per the TJ-analytic EFIT writer (handles the εa²·H₁ shift of the axis);
- evaluating ψ_plasma(r) from the radial ψ-ODE when r < 1, the TJ-analytic analytic vacuum solution (
GetPSIvacof Fitzpatrick's TJ) when 1 ≤ r < rc, and the 1/r² far-field form when r ≥ rc.
Reference: R. Fitzpatrick, TJ code (https://github.com/rfitzp/TJ) — the shape ODE (g₂, H₁, H₁', f₃), the GetPSIvac / GetHHvac vacuum extension, and the EFIT-writer (R, Z) → (r, w) Newton inversion that this routine adapts.
GeneralizedPerturbedEquilibrium.Equilibrium.tj_analytic_shape_initial — Method
Initial conditions at x = x0, matching the TJ-analytic model's near-axis expansion (cf. R. Fitzpatrick's TJ code, https://github.com/rfitzp/TJ).
GeneralizedPerturbedEquilibrium.Equilibrium.tj_analytic_shape_rhs! — Method
RHS for the TJ-analytic shape ODE (R. Fitzpatrick's TJ code parameterization, https://github.com/rfitzp/TJ). State: y[1]=ψ, y[2]=g₂, y[3]=H₁, y[4]=H₁', y[5]=f₃. The original derivation is written in x = r/a; we advance in physical r = a·x so d/dr = (1/a)·d/dx.
The params argument carries TJAnalyticShapeParams fields plus the current nu.
GeneralizedPerturbedEquilibrium.Equilibrium.tj_analytic_shape_solve — Method
Integrate the TJ-analytic shape ODE for the given ν. Pass saveat to collect output on a prescribed dense grid (used by tj_analytic_run_direct so the downstream Hₙ / ψ splines sit on uniform nodes); leave it nothing for the default adaptive save pattern used by tj_analytic_run.
Important types
EquilibriumConfig— top-level configuration container parsed from a
TOML file (outer constructor `EquilibriumConfig(path::String)` is
provided).EquilibriumControl— low-level control parameters (grid, jacobian
type, tolerances, etc.).PlasmaEquilibrium— the runtime structure containing spline fields,
geometry, profiles, and computed diagnostics (q-profile, separatrix,
etc.).LargeAspectRatioConfig,SolovevConfig— convenience structs to
construct analytic/model equilibria when using `eq_type = "lar"` or
`eq_type = "sol"`.Key functions
setup_equilibrium(path::String = "equil.toml")
— main entry point that reads configuration, builds the equilibrium,
runs the solver, and returns a `PlasmaEquilibrium` instance.equilibrium_separatrix_find!(pe::PlasmaEquilibrium)— locate
separatrix and related boundary geometry in-place.equilibrium_global_parameters!(pe::PlasmaEquilibrium)— populate
common global parameters (major radius, magnetic axis, volumes, etc.).equilibrium_qfind!(pe::PlasmaEquilibrium)— compute safety factor
(q) information across the grid.equilibrium_gse!(pe::PlasmaEquilibrium)— diagnostics on the
Grad–Shafranov solution.Example usage
Basic example: read a TOML config and build an equilibrium
using GeneralizedPerturbedEquilibrium
# Build from a TOML file (searches relative paths if needed)
pe = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium("docs/examples/ForceFreeStates.toml")
println("Magnetic axis: ", pe.params.r0, ", ", pe.params.z0)
println("q(0) = ", pe.params.q0)
# Find separatrix (in-place) and inspect results
GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_separatrix_find!(pe)
println("rsep = ", pe.params.rsep)Analytic / testing example: construct a large-aspect-ratio model
using GeneralizedPerturbedEquilibrium
# Create a LAR config from a small TOML fragment or file
larcfg = GeneralizedPerturbedEquilibrium.Equilibrium.LargeAspectRatioConfig(lar_r0=10.0, lar_a=1.0, beta0=1e-3)
pe = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(control=Dict("eq_filename"=>"unused","eq_type"=>"lar")), larcfg)
println("Built LAR equilibrium with a = ", larcfg.lar_a)Notes and Caveats
EquilibriumConfigis constructed from the[Equilibrium]section ofgpec.toml. Paths that are not absolute are resolved relative to the TOML file location.- The Equilibrium module contains readers for EFIT and CHEASE formats. Ensure the required data files are present and paths are set correctly in
gpec.toml.
See also
docs/src/stability.md— ideal MHD stability analysis built on top of the equilibrium