API Reference

solve(
    shot::Shot,
    its::Integer;
    tol::Real=0.0,
    relax::Real=1.0,
    debug::Bool=false,
    fit_fallback::Bool=true,
    concentric_first::Bool=true,
    concentric_last::Symbol=:warn,
    P::Union{Nothing,Tuple{<:Union{FE_rep,Function},Symbol},Profile}=nothing,
    dP_dψ::Union{Nothing,Tuple{<:Union{FE_rep,Function},Symbol},Profile}=nothing,
    F_dF_dψ::Union{Nothing,Tuple{<:Union{FE_rep,Function},Symbol},Profile}=nothing,
    Jt_R::Union{Nothing,Tuple{<:Union{FE_rep,Function},Symbol},Profile}=nothing,
    Jt::Union{Nothing,Tuple{<:Union{FE_rep,Function},Symbol},Profile}=nothing,
    Pbnd=shot.Pbnd,
    Fbnd=shot.Fbnd,
    Ip_target=shot.Ip_target
)

Solve the equilibrium, initially defined by shot with its iterations. Pressure and current information taken from shot unless provided in keywords

Returns a new Shot, often called refill by convention

Keyword arguments

• tol - Relative tolerance for convergence of the magnetic axis flux value to terminate iterations early
• relax - Relaxation parameter on the Picard iterations. Ψₙ₊₁ = relax * Ψ̃ₙ₊₁ + (1 - relax) * Ψₙ
• debug=true - Print debugging and convergence information
• fit_fallback=true - Use concentric surfaces if any flux surface errors on refitting. Improves robustness in early iterations
• concentric_first=true - Use concentric surfaces for first iteration, which can improve robustness if large changes from shot is expected
• concentric_last=:warn - If concentric surfaces used in final iteration, warn or error. Options are :warn, :error

FSA(f::F1, shot::F2, ρ::Real; tid=Threads.threadid()) where {F1,F2<:Shot}

Compute flux-surface average of f at ρ for the equilibrium defined in shot

Here f is a function of θ only, so something like f = θ -> F(ρ, θ) may be required

FSA(f::F1, shot::F2, ρ::Real, Vprime::F3; tid=Threads.threadid()) where {F1,F2<:Shot,F3<:FE_rep}

Compute flux-surface average of f at ρ for the equilibrium defined in shot, using the given finite-element representation of Vprime

Here f is a function of θ only, so something like f = θ -> F(ρ, θ) may be required

FSA(f::F1, shot::F2, ρ::Real, Vprime::Real; tid=Threads.threadid()) where {F1,F2<:Shot}

Compute flux-surface average of f at ρ for the equilibrium defined in shot, using the given value of Vprime

Here f is a function of θ only, so something like f = θ -> F(ρ, θ) may be required

Ip(shot::F1; ε::Real=1e-6) where {F1<:Shot}

Returns the plasma current of shot

Shot{
    I1<:Integer,
    VR1<:AbstractVector{<:Real},
    MR1<:AbstractMatrix{<:Real},
    MR2<:AbstractMatrix{<:Real},
    PT1<:Union{Nothing,Profile},
    PT2<:Union{Nothing,Profile},
    PT3<:Union{Nothing,Profile},
    PT4<:Union{Nothing,Profile},
    PT5<:Union{Nothing,Profile},
    R1<:Real,
    R2<:Real,
    IP1<:Union{Nothing,Real},
    FE1<:FE_rep,
    VFE1<:AbstractVector{<:FE_rep},
    Q1<:QuadInfo,
    VDC1<:Vector{<:DiffCache},
    F1<:Factorization
} <: AbstractEquilibrium

The fundamental data structure for a TEQUILA equilibrium, storing grid, flux-surface, and equilibrium information, as well as preallocated work arrays

shot(R,Z) returns the flux at point (R,Z)

VEQBoundary(; a, R0, Z0=0.0, B0, ka, c0=0.0, c_offsets=Float64[], s_offsets=Float64[])

Boundary MXH parameters (edge values of the shape profiles) in veqpy's convention: R = R0 + a(h + ρ cos θb), Z = Z0 + a(v − ρ k sin θ), θb = θ + c0 + Σ cm cos(mθ) + sm sin(mθ). c_offsets/s_offsets m=1-started.

VEQSource(; heat_profile, current_profile, Ip=NaN, beta=NaN)

PF-route source: heat_profile (pressure-gradient-like) and current_profile (FF'-like) sampled on a uniform normalized-flux axis with top.sample_count points. μ0 scaling of the pressure-like input happens at materialization, mirroring veqpy.

VEQTopology(; h_count, v_count=0, kappa_count, c0_count=0, psin_count,
              F_count=0, c_counts=Int[], s_counts=Int[], Nr, Nt,
              route=:PF, coordinate=:psin, nodes=:uniform,
              ip_constraint=false, beta_constraint=false,
              sample_count=51, K_max=nothing)

Static solve topology: number of radial Chebyshev coefficients per active profile family, quadrature grid size, and source route. c_counts/s_counts are m=1-started (like KernelTopology.s_counts in veqpy; note veqpy's ccounts are c0-started instead — here c0 has its own `c0count`).

Vprime(shot::F1, ρ::Real; tid=Threads.threadid()) where {F1<:Shot}

Compute dV/dρ at ρ for the equilibrium defined in shot

find_axis(shot::Shot)

Returns the location and flux value of the magnetic axis, (Raxis, Zaxis, Ψaxis)

find_axis(Ψ, R0::Real, Z0::Real)

Returns the location and flux value of the magnetic axis, (Raxis, Zaxis, Ψaxis), with initial guess (R0, Z0)

fsa_invR(shot::F1, ρ; tid=Threads.threadid()) where {F1<:Shot}

Compute <R⁻¹> at ρ for the equilibrium defined in shot

fsa_invR2(shot::F1, ρ; tid=Threads.threadid()) where {F1<:Shot}

Compute <R⁻²> at ρ for the equilibrium defined in shot

psi_ρθ(shot::Shot, ρ::Real, θ::Real)

Return the flux from shot at the MXH (ρ, θ)

veq_solve(kern::VEQKernel; x0=zeros(kern.x_size), tol=1e-9, maxiter=50)

Modified Newton on the packed residual with ForwardDiff Jacobian: the LU-factorized Jacobian is reused across iterations while the residual keeps contracting and refreshed when contraction stalls. Returns (x, converged, resnorm, iterations).

veq_solve!(shot::Shot;
    h_count=3, kappa_count=6, psin_count=6, c_counts=zeros(Int, length(shot.cfe)),
    s_counts=[3; zeros(Int, length(shot.sfe) - 1)],
    Nr=16, Nt=16, sample_count=51,
    outer_its=20, outer_tol=1e-8, tol=1e-9, debug=false)

Run the VEQ direct solve for shot's boundary and pressure/current profiles and write the result back into shot. All Shot profile routes are supported (dPdψ or P; FdFdψ, Jt, or JtR), on :poloidal or :toroidal grids, with optional shot.Ip_target enforced by TEQUILA-style current rescaling.

The flux scale ψs = −Ψaxis is converged by an outer fixed point (α2 = ψs²), warm-starting each inner Newton solve. When the profile conversion depends on the equilibrium (P/Jt/Jt_R, :toroidal grids, or Ip constraint), the solution is written back into shot every outer iteration so conversions and Ip see the current geometry.

Counts arrays may include trailing zeros: those harmonics stay passive (boundary value scaled by ρ^K_m) but still shape the geometry.

Returns shot.