MHD stability

VAFT answers stability questions at two levels of fidelity, and they are meant to be used in that order.

1. Analytic and empirical criteria — vaft.formula.stability. Closed-form limits (Greenwald, Troyon-type $\beta_N$, ballooning, kink, sawtooth) evaluated on scalars or NumPy arrays. Cheap enough to run on every time slice of every shot, so it is what the database summary uses to screen an operational space.
2. Linear MHD codes — vaft.code.gpec. VAFT prepares inputs for and drives the GPEC suite (DCON, RDCON, STRIDE, GPEC) from a reconstructed equilibrium, then collects the NetCDF products. This is where an actual $\delta W$ or $\Delta’$ comes from.

flowchart LR
    A["Reconstructed equilibrium<br/>(GEQDSK / ODS)"] --> B["vaft.formula.stability<br/>analytic limits"]
    A --> C["vaft.code.gpec<br/>prepare + run"]
    C --> D["DCON<br/>ideal delta W"]
    C --> E["RDCON / STRIDE<br/>tearing index Delta'"]
    C --> F["GPEC<br/>3-D perturbed response"]
    D --> G["NetCDF outputs<br/>collect_gpec_suite_outputs"]
    E --> G
    F --> G
    B --> H["stability_history.xlsx<br/>shot/time database summary"]
    G --> H

The closed-form criteria are catalogued with their exact signatures and unit traps on the Physics formulas page. This page covers what to do with them and how to run the codes.

Which mode family each tool addresses (interchange and Mercier, ballooning, kink and peeling, tearing and its parity channels), and why their criteria are not one severity scale, is on MHD mode representations across geometries.

What each criterion assumes and proves — exact definition, reduced model, empirical boundary, heuristic or solver-derived — is inventoried on Reduced stability diagnostics.

Screening an equilibrium

The stability functions are pure — they never touch an ODS — so pull the global quantities out yourself and hand them over.

import vaft

ods = vaft.omas.sample_ods()

# The packaged sample stores q_95 but not beta_normal; derive it (and beta_pol, beta_tor, li_3)
# from the stored equilibrium. A database shot whose equilibrium stage ran already carries it.
if 'equilibrium.time_slice.0.global_quantities.beta_normal' not in ods:
    vaft.omas.update_equilibrium_global_quantities_beta_li(ods, time_slice=0)

eq     = ods['equilibrium.time_slice.0.global_quantities']
beta_N = eq['beta_normal']        # [%·m·T/MA]
q_95   = eq['q_95']

B = vaft.formula.boundaries
beta_check = B.evaluate_boundary(B.get_boundary("troyon"), beta_N)   # beta_N <= 2.2 mu0 1e6 ~ 2.76
q_check    = B.evaluate_boundary(B.get_boundary("low_q"), q_95)      # q_psi > 2 (Greenwald 1988)

evaluate_boundary returns a margin that is positive on the permitted side whichever way the limit points, together with the source (paper, equation, DOI) and the inputs that fall outside the range the limit was fitted on. The Troyon coefficient is Troyon et al. (1984), p. 214: $(\beta A){max} \approx 2.2\,I_N$, i.e. $\beta_N \le 2.2\,\mu_0\cdot10^6 \approx 2.76$ %·m·T/MA, the origin of the rounded 2.8. kink_stability_criterion and beta_stability_boundary multiply the limit by $q{95}$, which no source does; they are deprecated (#350) and keep their old numbers only for compatibility.

The remaining local criteria return a (margin, critical_value) pair, and the margin is literally value - critical. A positive margin therefore means the plasma sits above the boundary:

alpha        = vaft.formula.ballooning_alpha_from_p_B_R(p, B, R, q)   # alpha = -2 mu0 R q^2 (dp/dR) / B^2
s            = vaft.formula.shear_from_r_q(r, q)                   # s = (r/q) dq/dr
d_alpha, a_c = vaft.formula.ballooning_stability_criterion(alpha, s)   # alpha_crit = 0.6 s

unstable = d_alpha > 0        # alpha has crossed the ideal ballooning boundary

Density limit

n_G = vaft.formula.greenwald_density(I_p=0.1, a=0.4)   # I_p in MA -> n_G in 1e19 m^-3
f_G = vaft.formula.greenwald_fraction(n_e=1.5, n_G=n_G)

greenwald_density takes $I_p$ in MA and returns $n_G$ in $10^{19}\ \mathrm{m^{-3}}$, and $f_G$ is conventionally formed with the line-averaged electron density. Both arguments to greenwald_fraction must use the same density definition and units — it just divides.

The combined margin helper (deprecated), and its one trap

plasma_stability_margins is deprecated (#350): its beta margin uses the unsourced $0.028\,q_{95}$ limit. Use evaluate_boundary with "troyon", "low_q" and "greenwald" instead, which share one sign convention. The description below documents the old behaviour.

beta_margin, q_margin, density_margin = vaft.formula.plasma_stability_margins(
    beta_N, q_95, n_e, n_G
)

The three returned numbers are not the same kind of quantity:

Note also that plasma_stability_margins is built on beta_stability_boundary, so it inherits the $0.028\,q_{95}$ convention, while kink_stability_criterion uses $2.8\,q_{95}$ — the same Troyon-type relation expressed in two different unit conventions. Do not mix the two in one plot.

Operational-space traces

The equilibrium IDS already carries $\beta_N$ and $q_{95}$ per time slice, so the whole-discharge view needs no formula call:

ods = vaft.database.load_ods(39915, source="public")

vaft.omas.plot_equilibrium_time_beta_n(ods)
vaft.omas.plot_equilibrium_time_q95(ods)
vaft.omas.plot_equilibrium_time_beta_p(ods)
vaft.omas.plot_equilibrium_time_beta_t(ods)

vaft.omas.plot_equilibrium_overview_histories(ods) packs $I_p$, $V_{\rm loop}$ and $\beta_N$ against $H_\alpha$, $B_z$ and $R_{\rm major}$ into one 3×2 figure — the fastest way to see whether a beta excursion coincides with a disruption.

Internal inductance against edge q

Two literature pictures look alike on this plane and are different objects (#1422):

What $l_i$ measures: the peaking of the current profile, through the poloidal-field energy inside the plasma normalised to the edge field. It does not fix the profile shape; many $j(r)$ give one $l_i$. Edge $q$ and $l_i$ together bound where the current gradient sits relative to the low-order rational surfaces, which is why both limits are drawn on this plane. Both papers are for conventional aspect ratio. At VEST aspect ratio, $q_\psi$, $q_{95}$ and the cylinder $q(a)$ differ substantially, and so do the $l_i$ conventions. A Wesson line is drawn only on a $q_\psi$ axis with $l_i(3)$, and a Cheng line only on the cylinder quantities. Neither is drawn on $q_{95}$.

import vaft
vaft.diagram.li_qa(reference="wesson_1989")      # empirical JET boundaries
vaft.diagram.li_qa(reference="cheng_1987")       # theoretical stable domain
from vaft.plot.operational_space import operational_space_population
fig, ax = operational_space_population(table, "li_qa_wesson")   # columns edge_safety_factor, internal_inductance_li3

Both curves were digitized at 600 dpi from the published figures, with the axes calibrated on the printed ticks. For Cheng, the printed closed form $\max(l_i/2) = [1 + 2\ln(q(a)/q(0))]/4$ is recovered to 0.01. The older vaft.formula.empirical_li_qa() arrays are the Wesson lower boundary of Fig. 6 (not Fig. 5, which is the Hugill diagram) and agree with the audit to about 0.03. Prefer the registered entries.

Atlas population histograms

The Tier A stability atlas population (#1852) is drawn as three stacked histograms, each against its own criterion: DCON ideal $W_t$ for $n = 1$–$6$ (unstable by sign only, $W_t < 0$; no $|W_t|$ band), each slice’s largest RDCON/STRIDE $\Delta’$ over its rational surfaces at $n = 1, 2$ (unstable if $\Delta’_{\max} > 0$), and the local Mercier $\max D_I$ (unstable $> 0$) and ballooning $\min C_A$ (unstable $< 0$). All values are the mpsi 512 run. The QA flags are counted in the reader’s summary table and do not filter the histograms.

from vaft.process.mhd_stability import load_stability_atlas, stability_atlas_populations
from vaft.plot.stability_atlas import stability_atlas_population
populations = stability_atlas_populations(*load_stability_atlas("~/runs/campaign/atlas/stability"))
print(populations.summary)           # slices, QA and unstable counts per series
fig, axes = stability_atlas_population(populations)

notebooks/stability_atlas_population.ipynb re-runs it from $VAFT_ATLAS_DIR.

Local MHD criteria from one equilibrium

vaft.process.local_mhd_criteria evaluates, from one reconstructed equilibrium and with no stability code:

• Mercier D_I and GGJ D_R (DCON’s formulas, checked against DCON’s own output);
• the normal and geodesic field-line curvature and the local shear, along flux surfaces and in R–Z;
• per-surface s–α in the #1720 normalisation (#1902).

The renderers in vaft.plot.local_mhd_criteria draw each against its own criterion: D > 0 is unstable, and the s–α boundary is the circular model’s, labelled as such.

from vaft.process.local_mhd_criteria import mercier_ggj_from_ods, s_alpha_from_ods
from vaft.plot.local_mhd_criteria import local_mhd_criteria_profile_interchange, local_mhd_criteria_s_alpha
profile = mercier_ggj_from_ods(ods, 0)            # D_I, D_R, H per surface
fig, ax = local_mhd_criteria_profile_interchange(profile)
fig, ax = local_mhd_criteria_s_alpha(s_alpha_from_ods(ods, 0))

Running the GPEC suite

vaft.code.gpec drives four external programs from one reconstructed GEQDSK. VAFT does not ship them: it writes their namelists, lays out the run directories, invokes the executables and collects the products.

The installation root is read from the GPECHOME environment variable unless you set gpec_home (or executable_dir) explicitly.

Configure, prepare, run

from pathlib import Path
from vaft.code import GPECCaseInputs, GPECSuiteConfig, run_gpec_suite_case

inputs = GPECCaseInputs(
    shot=39915,
    time_ms=340,
    geqdsk=Path("/srv/vest.filedb/public/39915/efit/g039915.00340"),
    workdir=Path("/srv/vest.filedb/public/39915/linear_stability"),
)

config = GPECSuiteConfig(
    modules=("dcon", "rdcon"),   # default: dcon, rdcon, stride, gpec
    modes=(1, 2),                # toroidal mode numbers n
    psilow=1e-2,
    psihigh=0.994,
    run_mode="run_if_available",
)

result = run_gpec_suite_case(inputs, config)
print(result.ok)

GPECSuiteConfig also exposes the DCON edge controls dcon_sas_flag (default False), dcon_qhigh (20.2) and dcon_psiedge (1.0), plus a per-module timeout (1200 s), an env mapping merged into the subprocess environment, and a backend that launches each module (default: a local child process; see execution backends).

run_mode selects what happens when the executables are not installed:

3-D coil input for any machine

By default the ideal-GPEC stage copies the packaged VEST coil.in template. Passing IdealGPECOptions(coil_specs=[CoilInputSpec(name, currents_a), ...]) instead generates coil.in and stages one <machine>_<set>.dat per activated set (coil_num is always the number of specs). For machine="vest" (the default) the geometry and the GPEC direction words come from the packaged VEST configuration and vaft.machine_mapping.conventions.VEST_GPEC_COIL_DIRECTIONS. For any other machine three things are mandatory and are never inherited from the template: coil_config, a mapping of set name to CoilSet3D (build one from Cartesian loops with vaft.machine_mapping.coils_non_axisymmetric_geometry.coil_set_from_xyz_loops, or from an existing GPEC file with coil_set_from_dat), and the machine’s ip_direction / bt_direction ("positive" for counter-clockwise, "negative" for clockwise, viewed from above, as GPEC defines them). A set built in memory is written out with the canonical header ncoil nsec npts nw (GPEC_COIL_DAT_HEADER). The toroidal-mode content of a current pattern and the vacuum field of the filaments are available as plain-array kernels in vaft.process.coils_non_axisymmetric.

To stage the inputs without executing anything — useful when the run itself is dispatched by a scheduler — call prepare_gpec_suite_case(inputs, config) directly. run_gpec(inputs, config) is a compatibility entry point that forwards to run_gpec_suite_case.

Directory layout and outputs

A prepared case materialises one directory per module and mode:

<workdir>/<time_label>/<module>/nn=<n>/

so the DCON $n=1$ run of the case above lands in .../linear_stability/00340/dcon/nn=1/. Collect the products of a whole case with:

from vaft.code import collect_gpec_suite_outputs

outputs = collect_gpec_suite_outputs("/srv/vest.filedb/public/39915/linear_stability")
for path in outputs["dcon"]:
    print(path)

The result is a dict keyed by module name (dcon, rdcon, stride, gpec).

GPECSuiteResult.records is a tuple of GPECModuleRun, one per module/mode, each carrying status (prepared, completed, stable, failed, skipped), returncode, reason, logs and outputs. The ok property is true when status is completed or stable and returncode == 0, so inspect it per module rather than trusting the suite-level return code alone:

for record in result.records:
    print(record.module, record.mode, record.status, record.ok, record.reason)

stable is a successful status, kept apart from completed because a stable equilibrium produces no unstable-mode output and so cannot be told from a broken run by what it wrote (#423). It is read from the solver’s own free-boundary energy – Re(total1), which DCON, RDCON and STRIDE all write into their netCDF – so a companion such as rmatch exiting badly with nothing to analyse does not make the cell a failure. A run that produced no usable output is still failed, stable or not.

GEQDSK headers

GPEC’s EFIT reader is fixed-column and rejects headers that EFIT itself will happily write. If a case fails at the file-reading stage, normalise the header first:

from vaft.code import format_gfile_header_for_gpec

header = format_gfile_header_for_gpec(header_line, shot=39915, time_ms=340)

It recovers the shot and time from the existing header when you do not pass them.

Batch runs and the stability database

Running one case at a time is only useful for debugging. Two Snakemake workflows do the real work.

Sweeping shots and times — workflow/running_linear_stability/ discovers the reconstructed g-files under BASE_DIR (/srv/vest.filedb/public by default, set in config.yaml) and fans out the codes over the mode scans: $n = 1\ldots6$ for the ideal (DCON) scan and $n = 1, 2$ for the resistive (RDCON/STRIDE) scan.

cd workflow/running_linear_stability
snakemake --cores 16

Mining the results into a table — gen_stability_history.py walks {shot}/linear_stability/{time}/{dcon,rdcon,stride}/nn={n}/, reads the NetCDF products and emits stability_history.xlsx, one row per (shot, time):

The verdicts follow the standard linear-MHD reading, with a guard against numerical blow-ups:

• Ideal: "Stable" when $\delta W > 0$ for every scanned $n$; "Unstable" if any $n$ has $\delta W < 0$. If any $|\delta W| \geq 100$ the row is flagged "Error" rather than trusted.
• Resistive: "Stable" when no scanned $\Delta’$ is positive; "Unstable" as soon as one is. Any $|\Delta’| \geq 1000$ flags the row "Error".

Because a single bad eigenvalue poisons the verdict, always filter on ideal_stability != "Error" before drawing conclusions from the sheet.

Companion scripts in the same directory join the stability table to the CHEASE equilibrium history (join_chease_stability_history.py, merging on shot + time [ms]), plot the resulting history (plot_stability_history.py), and draw the operational-limit diagram (plot_stability_limit.py). See Signal processing and EM modeling for how these pipelines fit together with the rest of the automatic processing.

Notebooks

Three notebooks in the repository frame the stability workflow. They are currently outlines — they define the inputs, configuration conventions and expected products for each code rather than executing them, so read them as specifications, not as runnable examples.

• linear_ideal_stability_analysis_with_dcon.ipynb
• linear_resistive_stability_analysis_with_rdcon.ipynb
• perturbed_equilibrium_and_3d_response_with_gpec.ipynb

For equilibria that actually run end to end today, start from mhd_equilibrium_analysis.ipynb and equilibrium_refinement_using_chease.ipynb, which produce the refined equilibria the GPEC suite consumes — see the Examples page.

See also

• Physics formulas — full signature reference for vaft.formula.stability, including the beta conversions, characteristic speeds and unit traps.
• Equilibrium — producing the equilibrium the codes consume.
• Ballooning formulations — how the reduced $s$–$\alpha$ model, DCON’s $C_A$ and GPEC.jl’s ballooning $\Delta’$ relate, and the shared normalisation.
• MHD mode representations across geometries — how slab, cylindrical and toroidal mode families relate, and which tool sits where.
• Reduced stability diagnostics — the logical status, assumptions and validation path of every analytic and reduced criterion.
• Data structures (ODS, IDS, IMAS) — where beta_normal, q_95 and the rest of equilibrium.time_slice[:].global_quantities live.
• Source: vaft/formula/stability.py and vaft/code/gpec.py.