GitHub

Scripting API

GPEC's pipeline is normally driven from a gpec.toml deck through GeneralizedPerturbedEquilibrium.main. The scripting API exposes the same stages as ordinary Julia functions, so a run can be built, parameterized and looped over in a script without writing a deck.

using GeneralizedPerturbedEquilibrium

eq  = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada")
ffs = solve(eq, Riccati(); nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true)
rmp = RMPField("coils.dat")
pe  = perturbed_equilibrium(ffs, rmp)

solve(eq, alg; kwargs...) is sugar for the canonical problem form: a EulerLagrangeProblem names WHAT is being solved (the perturbed-plasma Euler-Lagrange system posed on this equilibrium, with this mode range, wall and closure) and the integrator names HOW. A PlasmaEquilibrium hosts many possible problems; the problem type keeps solve unambiguous as other problem classes appear.

prob = EulerLagrangeProblem(eq; nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true)
ffs  = solve(prob, Riccati())

The four objects map one-to-one onto the pipeline stages:

  • PlasmaEquilibrium reads and processes the equilibrium — the [Equilibrium] section. Analytic equilibria (sol, lar, tj_analytic) take their parameters from a separate config object and are built with setup_equilibrium(eq_config, analytic_config) instead.
  • solve runs the force-free-states stage — the [ForceFreeStates] section — and returns a ForceFreeStatesResult, exactly the object the TOML driver publishes.
  • RMPField describes the external field — the [ForcingTerms] section — without reading anything from disk until it is applied.
  • perturbed_equilibrium runs the plasma-response stage — the [PerturbedEquilibrium] section.

Combining forcing sources

RMPFields form a vector space: +, - and scalar * build lazy linear combinations that record their terms and compute nothing until the perturbed-equilibrium stage materializes them. Because the plasma response is linear in the forcing, driving with a combination is exactly equivalent to combining the individually materialized fields — each term is evaluated on the control surface and the mode amplitudes are summed. A complex scalar phase-rotates a source.

nominal    = RMPField("nominal_efc.dat")
weld_field = RMPField("weld_fields.h5")
pe = perturbed_equilibrium(ffs, 2.0 * nominal + 0.5 * weld_field)

This is the substrate for error-field workflows: build per-unit sources independently (a coil set per ampere, a displacement field per millimeter), then assemble physical cases by weighting and summing — without recomputing anything per combination until the final perturbed_equilibrium call.

The weights are linear-combination coefficients, not physical amplitudes: sources that are not scalar multiples of each other get their own description and the algebra. A coil set with a failed conductor is nominal - failed_coil, not 0.9 * nominal.

Choosing an integrator

The second argument of solve picks the formalism, and its fields are that formalism's tunables. Everything else is a ForceFreeStatesControl keyword, so the TOML keys and the solve keywords are the same knobs:

ffs = solve(eq, Forward(); nn=1, vac_flag=true)                 # dense ξ profiles for PerturbedEquilibrium
ffs = solve(eq, Riccati(; nchunks=40); nn=1, vac_flag=true)     # chunked propagators, Δ′ matrix
ffs = solve(eq, Galerkin(; nx=512); nn=1)                       # RDCON outer-region Galerkin Δ′

Which products each formalism can supply differs; a result carries nothing in the fields its integrator does not produce and consumers warn and skip rather than erroring. See the Stability Analysis page for the result struct and its capability gates.

Inner-layer matching is requested with the integrator-agnostic match keyword, which closes the basis with a resistive layer solution instead of the ideal jump condition:

ffs = solve(eq, Galerkin(); nn=1,
    match=ResistiveMatch(; eta=[1e-6, 2e-6], rho=[1e-7, 1e-7], rotation=[0.0, 0.0]))
@assert ffs.closure === :matched

Only the Galerkin formalism implements the match today; requesting one from Forward or Riccati errors. Kinetic runs (kinetic_factor > 0) need the [KineticForces] profiles and remain TOML-driven.

Entry points

GeneralizedPerturbedEquilibrium.EulerLagrangeProblem — Type
EulerLagrangeProblem(equil; nn, wall=Vacuum.WallShapeSettings(), match=nothing,
                     dir_path=".", debug=DebugSettings(), kwargs...)

The perturbed-plasma Euler-Lagrange problem posed on an equilibrium: the extremization of the perturbed potential energy whose solutions are the force-free (and, via the TOML path, kinetic) perturbed states. This is the WHAT of a stability solve; the integrator passed to solve is the HOW. A PlasmaEquilibrium hosts many possible problems — this type names this one, so solve stays unambiguous as other problem classes appear.

nn is the toroidal mode number or range. wall is the vacuum wall shape, match an optional ResistiveMatch closing the basis with an inner-layer solution instead of the ideal jump, dir_path the working directory outputs are written to, and debug the diagnostic dump settings of the DEBUG deck section. Any remaining keyword is a ForceFreeStatesControl field, so the TOML keys and the problem keywords are the same knobs. nn_low/nn_high are rejected — they come from nn.

Fields

  • equil::Equilibrium.PlasmaEquilibrium - The equilibrium the problem is posed on.
  • wall::Vacuum.WallShapeSettings - Vacuum wall shape for the free-boundary energies.
  • match::Union{Nothing,ForceFreeStates.ResistiveMatch} - Optional inner-layer closure.
  • dir_path::String - Working directory for outputs.
  • debug::DebugSettings - Diagnostic dump settings.
  • ctrl_kwargs::Dict{Symbol,Any} - ForceFreeStatesControl keywords, nn already folded in.
source
CommonSolve.solve — Function
solve(prob::EulerLagrangeProblem, alg) -> ForceFreeStatesResult
solve(equil, alg; nn, kwargs...) -> ForceFreeStatesResult

Solve the perturbed-plasma EulerLagrangeProblem with the formalism alg — Forward, Riccati or Galerkin — and return the published ForceFreeStatesResult. This is the scripting entry point; it runs the same stages a gpec.toml run of main does and produces the same result object. The second form is sugar building the problem from an equilibrium and the problem keywords in one call.

Knobs owned by alg or match are rejected as ForceFreeStatesControl keywords. Kinetic runs are TOML-driven this cycle: kinetic_factor > 0 needs the [KineticForces] profiles and errors here.

eq   = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada")
prob = EulerLagrangeProblem(eq; nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true)
ffs  = solve(prob, Riccati())
ffs  = solve(eq, Riccati(); nn=1, vac_flag=true)   # equivalent one-line form
source
solve(equil::PlasmaEquilibrium, alg; nn, kwargs...) -> ForceFreeStatesResult

Convenience form of solve: builds the EulerLagrangeProblem from the equilibrium and the problem keywords, then solves it with alg.

source
GeneralizedPerturbedEquilibrium.perturbed_equilibrium — Function
perturbed_equilibrium(ffs, rmp; forcing_modes=nothing, coil_sets=nothing, kwargs...) -> PerturbedEquilibriumState

Compute the plasma response to the external field rmp on top of a force-free-states solve ffs, and write the perturbed-equilibrium outputs. rmp is an RMPField; keyword arguments are PerturbedEquilibrium.PerturbedEquilibriumControl fields.

Products the producing integrator could not supply gate the corresponding calculation: the step warns and is skipped rather than erroring, so a Riccati- or Galerkin-fed result still flows through.

forcing_modes injects already-loaded modes (the gpec.h5 replay path) and coil_sets already-built coil geometry, both bypassing the corresponding read.

pe = perturbed_equilibrium(ffs, RMPField("forcing.dat"))
source

Integrator selectors

GeneralizedPerturbedEquilibrium.ForceFreeStates.AbstractIntegrator — Type
AbstractIntegrator

Supertype of the three force-free-states formalisms selected by the scripting API solve(equil, alg; ...): Forward, Riccati and Galerkin.

An integrator object is pure configuration. solve translates it into the matching ForceFreeStatesControl keywords, so the struct fields and the TOML [ForceFreeStates] keys always describe the same solve — ForceFreeStatesControl stays the single source of truth and the TOML path is unaffected.

source
GeneralizedPerturbedEquilibrium.ForceFreeStates.Forward — Type
Forward()

Serial Euler-Lagrange integrator: sweeps the full radial domain and stores the dense ξ solution. The only formalism that supports kinetic runs and the only one whose solution feeds the profile-based PerturbedEquilibrium outputs. Maps onto integrator = "forward".

source
GeneralizedPerturbedEquilibrium.ForceFreeStates.Galerkin — Type
Galerkin(; solver="LU", nx=256, ...)

RDCON outer-region singular Galerkin solver: solves the same Euler-Lagrange system variationally on a finite-element grid packed around the rational surfaces, producing Δ′ without a radial ODE sweep. Maps onto integrator = "galerkin"; every field is the matching gal_* control key without the prefix.

Fields

  • solver::String - Banded linear solver, "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). Matching requires "LU".
  • nx::Int - Elements per interval between singular surfaces.
  • nq::Int - Gauss-Lobatto quadrature order per element.
  • pfac::Float64 - Grid packing ratio near singular surfaces.
  • dx0::Float64 - Resonant-element integration truncation distance, in units of 1/|n q′|.
  • dx1::Float64 - Resonant-element size, in units of 1/|n q′|.
  • dx2::Float64 - Extension-element size, in units of 1/|n q′|.
  • cutoff::Int - Number of elements carrying the large solution as driving term.
  • tol::Float64 - Resonant-quadrature tolerance.
  • gnstep::Int - Maximum resonant-quadrature evaluations.
  • dx1dx2_flag::Bool - Enable the special dx1/dx2 treatment of resonant and extension elements.
  • sing_order::Int - Base power-series order for the singular asymptotics.
  • sing_order_ceiling::Bool - Auto-raise the order per surface for a high Mercier index.
  • rpec_flag::Bool - Append the mpert coil-response columns to the Δ′ solve. Forced on when a ResistiveMatch is requested.
  • edge_onesided::Bool - Pack the two end intervals one-sided toward their single rational end instead of the Fortran symmetric pack.
source
GeneralizedPerturbedEquilibrium.ForceFreeStates.ResistiveMatch — Type
ResistiveMatch(; eta=[], rho=[], rotation=[], gamma=5/3, ideal=false, inner_solver="ray", ...)

Inner-layer matching configuration, passed to solve as match= and independent of the integrator that produced the outer solution. Requesting a match closes the basis with a resistive inner-layer solution instead of the ideal jump condition, so the result carries closure = :matched and a non-zero bpen.

Only Galerkin implements the match today; a Riccati or Forward solve with match set errors. The per-surface vectors are ordered core to edge and must have one entry per matched rational surface.

Fields

  • eta::Vector{Float64} - Per-surface resistivity η.
  • rho::Vector{Float64} - Per-surface mass density ρ in kg/m³.
  • rotation::Vector{Float64} - Per-surface rotation frequency f in Hz; the forced eigenvalue is γ_s = 2πi·n·f.
  • gamma::Float64 - Ratio of specific heats Γ in the resistive-layer coefficients.
  • ideal::Bool - Build the ideal (perfectly shielded) matched solution: skip the inner layer and use the bare coil columns. eta, rho and rotation are then unread.
  • inner_solver::String - Inner-layer Δ backend, "ray" (rotated-contour collocation) or "galerkin" (Hermite-cubic elements).
  • inner_xfac::Float64 - Asymptotic-matching radius multiplier of the "galerkin" backend.
  • inner_nx::Int - Grid cells of the "galerkin" backend.
  • inner_nq::Int - Quadrature order per cell of the "galerkin" backend.
  • inner_cutoff::Int - Cells carrying the large solution as driving term in the "galerkin" backend.
  • inner_kmax::Int - Large-x asymptotic series order of the "galerkin" backend.
source
GeneralizedPerturbedEquilibrium.ForceFreeStates.Riccati — Type
Riccati(; nchunks=0)

STRIDE-style chunked Riccati integrator: solves the Euler-Lagrange system on independent radial chunks and couples them through a boundary-value problem, which is what unlocks the inter-surface Δ′ matrix. Threads come from julia -t; the chunk count is the only tunable and never depends on the thread count. Maps onto integrator = "riccati".

Fields

  • nchunks::Int - Chunk-count target; 0 derives it from the number of singular surfaces.
source