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_typeMethodBest for
"efit"Geometric-angle field-line ODEDefault; fast and robust for standard cases
"efit_arclength"Arc-length field-line ODENear-separatrix surfaces where the geometric-angle ODE becomes singular
"efit_by_inversion"Contour.jl marching-squares → inverse solverHighest 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.AnalyticEqSpecType
AnalyticEqSpec

Single-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::Stringgpec.toml section key (e.g. "SOL_INPUT")
  • config_type::DataType*Config type, constructed via its (::Dict) / (::String) ctor
  • run_fn::Function(EquilibriumConfig, config) -> eq_input solver entry point
source
GeneralizedPerturbedEquilibrium.Equilibrium.DirectIngestType
DirectIngest

The 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 spline
  • sq_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)
source
GeneralizedPerturbedEquilibrium.Equilibrium.DirectRunInputType
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::EquilibriumConfig The equilibrium configuration object.

  • sq_in 1D spline data versus normalized poloidal flux psin. Quantities:

    1. F = R * B_t — toroidal flux function [m·T]
    2. μ₀ * Pressure — plasma pressure (non-negative) [T²]
    3. q — safety factor profile
    4. √ψ_norm — square root of normalized flux
  • psi_in 2D 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:

    1. ψ(R, Z) = ψ_boundary - ψ(R, Z)

    2. ψ = ψ * sign(ψ(centerR, centerZ))

      • 1D profiles are represented by CubicInterpolant or CubicSeriesInterpolant
      • 2D flux surfaces by CubicInterpolantND
  • 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-files

  • ingest::EquilibriumIngest — captured raw arrays for the gpec.h5 rerun snapshot (a DirectIngest for file-based reads, or nothing for analytic equilibria)

source
GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfigType
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 file
  • jac_type::String - Jacobian coordinate type ("hamada", "pest", "equal_arc", "boozer", "park", "other")
  • power_bp::Int - Poloidal field power exponent for Jacobian
  • power_b::Int - Total field power exponent for Jacobian
  • power_r::Int - Major radius power exponent for Jacobian
  • power_rc::Int - Minor radius (rfac = √((R-R₀)²+(Z-Z₀)²)) power exponent for Jacobian
  • r0exp::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 coordinate
  • psihigh::Float64 - Upper limit of normalized flux coordinate
  • mpsi::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 points
  • newq0::Int - Override for on-axis safety factor (0 = use input value)
  • etol::Float64 - Error tolerance for equilibrium solver
  • force_termination::Bool - Terminate after equilibrium setup (skip stability calculations)
  • use_galgrid::Bool - Use the same grid as galerkin method
source
GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumParametersType
EquilibriumParameters

A 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 location
  • b0::Union{Nothing,Float64} - Total magnetic field strength at the axis [T]
  • q0::Union{Nothing,Float64} - Safety factor at the axis
  • qmin::Union{Nothing,Float64} - Minimum safety factor in the plasma
  • qmax::Union{Nothing,Float64} - Maximum safety factor in the plasma
  • qa::Union{Nothing,Float64} - Safety factor at the plasma edge
  • q95::Union{Nothing,Float64} - Safety factor at 95% flux surface
  • qextrema_psi::Union{Nothing,Vector{Float64}} - Normalized flux values at q extrema
  • qextrema_q::Union{Nothing,Vector{Float64}} - Safety factor values at extrema
  • mextrema::Union{Nothing,Int} - Number of extrema in q-profile
  • psi_norm::Union{Nothing,Float64} - Normalized poloidal flux
  • b_norm::Union{Nothing,Float64} - Normalized magnetic field strength
  • psi_axis::Union{Nothing,Float64} - Poloidal flux at the axis
  • psi_boundary::Union{Nothing,Float64} - Poloidal flux at the boundary
  • psi_boundary_norm::Union{Nothing,Float64} - Normalized boundary flux
  • psi_axis_norm::Union{Nothing,Float64} - Normalized axis flux
  • psi_boundary_offset::Union{Nothing,Float64} - Boundary flux offset
  • psi_axis_offset::Union{Nothing,Float64} - Axis flux offset
  • psi_boundary_sign::Union{Nothing,Int} - Sign of boundary flux
  • psi_axis_sign::Union{Nothing,Int} - Sign of axis flux
  • psi_boundary_zero::Union{Nothing,Bool} - Whether boundary flux is zero
  • rmean::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 elongation
  • delta1::Union{Nothing,Float64} - Upper triangularity
  • delta2::Union{Nothing,Float64} - Lower triangularity
  • bt0::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 output
  • diagnose_src::Bool - Enable source data diagnostics
  • diagnose_maxima::Bool - Enable extrema diagnostics
  • volume::Union{Nothing,Float64} - Plasma volume [m³]
  • betat::Union{Nothing,Float64} - Toroidal beta
  • betan::Union{Nothing,Float64} - Normalized beta
  • betaj::Union{Nothing,Float64} - Total beta
  • betap1::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)
source
GeneralizedPerturbedEquilibrium.Equilibrium.GeometryProfileSplinesType
GeometryProfileSplines

Named 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 (matches rzphi_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]
source
GeneralizedPerturbedEquilibrium.Equilibrium.InverseIngestType
InverseIngest

The 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 spline
  • sq_fs::Matrix{Float64} — 1D profile node values
  • rz_xs::Vector{Float64} — ψ grid for the R, Z maps
  • rz_ys::Vector{Float64} — θ grid for the R, Z maps
  • R_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]
source
GeneralizedPerturbedEquilibrium.Equilibrium.InverseRunInputType
InverseRunInput(...)

A container struct for inputs to the inverse_run function.

Fields

  • config::EquilibriumConfig - The equilibrium configuration object
  • sq_in::CubicSeriesInterpolant - 1D profile spline for F, P, q
  • rz_in_xs::Vector{Float64} - ψ coordinate grid for rz_in
  • rz_in_ys::Vector{Float64} - θ coordinate grid for rz_in
  • rz_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 the gpec.h5 rerun snapshot (an InverseIngest for file-based reads, or nothing for analytic equilibria)
source
GeneralizedPerturbedEquilibrium.Equilibrium.KineticProfileSplinesType
KineticProfileSplines

Named 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 logarithm
  • nui_spline, nue_spline: Krook collision frequencies [s⁻¹]
  • zeff_spline: Effective charge Z_eff
  • ni_deriv, ne_deriv, Ti_deriv, Te_deriv: Derivative views for ψ-derivative access
source
GeneralizedPerturbedEquilibrium.Equilibrium.KineticProfileSplinesMethod
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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.LargeAspectRatioConfigType
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 as p00 * (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 as sigma0 * (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
source
GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibriumType
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 via profiles.F_spline.y[i], etc. Access derivatives via profiles.F_deriv.y[i] or profiles.F_deriv(psi).

  • geometry::GeometryProfileSplines: Named 1D splines for flux-surface-averaged geometry (area, ⟨r⟩, ⟨R⟩), populated automatically by compute_geometry_profiles during 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 the gpec.h5 rerun snapshot — a DirectIngest/InverseIngest for file-based equilibria, or nothing for analytic ones (regenerated from their TOML section on replay)

source
GeneralizedPerturbedEquilibrium.Equilibrium.ProfileSplinesType
ProfileSplines

Named 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 xs field or spline.cache.x
source
GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfigType
SolovevConfig(...)

A mutable struct holding parameters for the Solev'ev (SOL) plasma equilibrium model.

Fields:

  • mr: number of radial grid zones
  • mz: number of axial grid zones
  • ma: number of flux grid zones
  • e: elongation
  • a: minor radius
  • r0: major radius
  • q0: safety factor at the o-point
  • p0fac: 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)
source
GeneralizedPerturbedEquilibrium.Equilibrium.TJAnalyticConfigType
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

source
GeneralizedPerturbedEquilibrium.Equilibrium.TJAnalyticShapeParamsType

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
source
GeneralizedPerturbedEquilibrium.Equilibrium._build_psi_gridMethod
_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.

source
GeneralizedPerturbedEquilibrium.Equilibrium._cubic_resampleMethod

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

source
GeneralizedPerturbedEquilibrium.Equilibrium._estimate_log_slopeMethod
_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.

source
GeneralizedPerturbedEquilibrium.Equilibrium._estimate_mid_spacingMethod
_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].

source
GeneralizedPerturbedEquilibrium.Equilibrium._read_1d_gfile_formatMethod

read1dgfileformat(linesblock, numvalues)

Internal helper function to parse Fortran-style fixed-width numerical blocks from a vector of strings.

Arguments:

  • lines_block: A Vector{String} containing the lines to parse.
  • num_values: The total number of Float64 values to read from the block.

Returns:

  • A Vector{Float64} containing the parsed values.
source
GeneralizedPerturbedEquilibrium.Equilibrium._read_kinetic_tableMethod

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.

source
GeneralizedPerturbedEquilibrium.Equilibrium._winding_number_vertsMethod
_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).

source
GeneralizedPerturbedEquilibrium.Equilibrium.adaptive_grid_paramsMethod
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(β).

source
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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.arclength_fieldline_intMethod
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/Bp
  • y_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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.build_direct_from_ingestMethod
build_direct_from_ingest(config::EquilibriumConfig, ingest::DirectIngest) -> DirectRunInput

Rebuild 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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.classify_topologyMethod
classify_topology(raw_profile, psio; xpt_threshold=0.05) → Symbol

Classify 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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.compute_geometry_profilesMethod
compute_geometry_profiles(rzphi_xs, rzphi_ys,
                          rzphi_rsquared, rzphi_offset, rzphi_jac, ro)
    → GeometryProfileSplines

Build 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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.compute_sqrtamatMethod
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.

source
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: A FieldLineDerivParams struct with all necessary parameters.
  • eta: The independent variable (geometric angle η).
source
GeneralizedPerturbedEquilibrium.Equilibrium.direct_fieldline_intMethod
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: DirectRunInput object 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: A DirectBField object with values at the integration start point.
source
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 mutable DirectBField struct to store the results
  • r: R-coordinate to evaluate at
  • z: Z-coordinate to evaluate at
  • psi_in: 2D cubic interpolant for poloidal flux ψ(R,Z)
  • sq_in: 1D cubic spline for profiles F(ψ_norm) and P(ψ_norm)
  • sq_in_deriv: Pre-computed derivative view of sq_in
  • psio: total toroidal flux
  • derivs: An integer specifying number of derivatives to compute (0, 1, or 2)
source
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: A DirectRunInput object 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)
source
GeneralizedPerturbedEquilibrium.Equilibrium.direct_refineMethod
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: A FieldLineDerivParams struct.

Returns:

  • The refined rfac value.
source
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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_solverFunction
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: A DirectRunInput object containing the initial splines (psi_in, sq_in) and run parameters (equil_input).

Returns:

  • A PlasmaEquilibrium object containing the final, processed equilibrium data, including the profile spline (sq), the coordinate mapping spline (rzphi), and the physics quantity spline (eqfun).
source
GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_solver_by_inversionMethod
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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.flux_surface_metricMethod
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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.inverse_extrapMethod
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 corresponding xx.
  • x: A scalar Float64 value at which to extrapolate.

Returns:

  • A vector of length n representing the extrapolated function values at x.
source
GeneralizedPerturbedEquilibrium.Equilibrium.lar_init_conditionsMethod
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 of lar_a).
  • lar_input: A LargeAspectRatioConfig object containing equilibrium parameters.

Returns:

  • r: Physical radius corresponding to rmin * lar_a.
  • y: Initial state vector of length 5.
source
GeneralizedPerturbedEquilibrium.Equilibrium.load_kinetic_profilesMethod
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)
    → KineticProfileSplines

Parse 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 file
  • zi, zimp: Main ion and impurity charge numbers
  • mi, 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 normalization 2π·ψ₀ — required when density_factor, temperature_factor, or toroidal_rotation_factor differ from 1.0

Scaling sequence

When any of density_factor, temperature_factor, toroidal_rotation_factor differ from 1.0:

  1. Build first-pass cubic splines from unscaled profiles (for derivatives)
  2. Compute diamagnetic frequencies wdian, wdiat and total toroidal rotation wphi
  3. Scale: wdian_new = temperature_factor * wdian, wdiat_new = temperature_factor * wdiat (density_factor cancels in T*(dn/dψ)/n)
  4. Reform: omegaE = toroidal_rotation_factor * wphi - wdian_new - wdiat_new
  5. 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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.make_multi_stretched_gridMethod
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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.make_optimal_mpsiMethod
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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.make_optimal_psi_gridMethod
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 ψ.

source
GeneralizedPerturbedEquilibrium.Equilibrium.make_stretched_z_gridMethod
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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.read_chease_asciiMethod
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: The EquilibriumConfig object containing the filename and parameters.

Returns:

  • A InverseRunInput object ready for the inverse solver.
source
GeneralizedPerturbedEquilibrium.Equilibrium.read_efitMethod
_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: The EquilInput object containing the filename and parameters.

Returns:

  • A DirectRunInput object ready for the direct solver.
source
GeneralizedPerturbedEquilibrium.Equilibrium.read_imasMethod
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: EquilibriumConfig with eq_type = "imas" and imas_cocos set.
  • dd: populated IMASdd.dd with dd.equilibrium.time_slice[] containing:
    • global_quantities.psi_axis, global_quantities.psi_boundary
    • profiles_1d.psi, profiles_1d.f, profiles_1d.pressure, profiles_1d.q
    • profiles_2d[1].grid.dim1 (R), profiles_2d[1].grid.dim2 (Z), profiles_2d[1].psi
source
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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.rootarea_to_area_weightMethod
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 to the area-weighted field : 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 () basis is the one in which operator singular values / spectra are independent of the straight-field-line (working) coordinate — see field_space_response_matrices. is not coordinate-invariant; it is only a field/flux recovery view. [Pharr 2026]

source
GeneralizedPerturbedEquilibrium.Equilibrium.sol_runMethod

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 zones
  • mz: Number of axial grid zones
  • ma: Number of flux grid zones
  • e: Elongation
  • a: Minor radius
  • r0: Major radius
  • q0: Safety factor at the o-point
  • p0fac: 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:

  • DirectRunInput object
source
GeneralizedPerturbedEquilibrium.Equilibrium.tj_analytic_f1Method
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

source
GeneralizedPerturbedEquilibrium.Equilibrium.tj_analytic_find_nuMethod

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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.tj_analytic_runMethod
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/a

with ν = 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

source
GeneralizedPerturbedEquilibrium.Equilibrium.tj_analytic_run_directMethod
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 (GetPSIvac of 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.

source
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.

source
GeneralizedPerturbedEquilibrium.Equilibrium.tj_analytic_shape_solveMethod

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.

source

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

  • EquilibriumConfig is constructed from the [Equilibrium] section of gpec.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