spectre.converters.vmec_helpers

VMEC physics helpers and generic VMEC-to-SPECTRE builder.

Contains everything that is shared across VMEC adapters:

Adding a new VMEC interface

Write a new adapter module that imports ContinuousPhysics from here and implements its own extraction function:

def extract_<name>_physics(source) -> ContinuousPhysics:
    ...
    return ContinuousPhysics(
        curtor=..., phiedge=..., nfp=...,
        rbc={...}, zbs={...},
        profile_type='iota',   # or 'current'
        pressure=...,          # callable pressure(s), mu0-scaled
        iota=...,              # callable iota(s)
        current_density=...,   # callable j(s) = dI/ds
    )

Then call build_spectre_input() directly — no other file in this package needs to change.

Module Contents

spectre.converters.vmec_helpers.MU0
spectre.converters.vmec_helpers.power_series(coeffs)

Return a callable that evaluates sum_i coeffs[i] * s**i.

VMEC stores pressure, iota, and current profiles as power-series coefficients in normalised toroidal flux s ∈ [0, 1].

Args:

coeffs: sequence of polynomial coefficients ordered from degree 0.

Returns:

Callable f(s) that accepts scalars or arrays.

spectre.converters.vmec_helpers.pressure_function(am, pres_scale)

Return a callable pressure(s) for VMEC pressure coefficients.

The returned callable evaluates the am power series scaled by pres_scale and mu0, so the result is in SPECTRE’s internal pressure units.

Args:

am: VMEC pressure polynomial coefficients (array-like). pres_scale: VMEC pressure scale factor (multiplies the polynomial).

Returns:

Callable pressure(s) accepting scalars or arrays of normalised toroidal flux s ∈ [0, 1].

spectre.converters.vmec_helpers.compute_pressure_profile(pressure, tfluxes)

Volume-averaged pressure in each sub-volume between consecutive interfaces.

SPECTRE stores a single pressure value per volume. We obtain it by integrating the continuous profile pressure(s) over the s-interval [s_{i-1}, s_i] of each volume and dividing by the interval width, giving the volume average.

Args:
pressure: callable pressure(s) returning the pressure in SPECTRE

units (i.e. already mu0-scaled) at normalised toroidal flux s.

tfluxes: 1-D array of interface toroidal-flux values in (0, 1].

Returns:

1-D array of length len(tfluxes) with <pressure> per volume.

spectre.converters.vmec_helpers.iota_function(vmec)

Return a callable iota(s) for a simsopt Vmec object.

Iota is the one profile whose representation differs by backend, so the dispatch lives here:

  • indata source (loaded from input.*): evaluate the ai power-series coefficients via power_series(), matching VMEC’s own profile evaluation.

  • wout source (loaded from wout_*.nc): cubic interpolation of the converged full-mesh iotaf grid. The ai polynomial fit may not reproduce iotaf exactly, so the output grid is preferred when present.

Args:
vmec: simsopt Vmec instance with either an indata (input file)

or wout (netCDF output) attribute.

Returns:

Callable iota(s) accepting scalars or arrays of normalised toroidal flux s ∈ [0, 1].

Raises:

RuntimeError: if vmec has neither indata nor wout.

spectre.converters.vmec_helpers.normalized_enclosed_current(current_density, s)

Normalised cumulative enclosed toroidal current Î(s) at flux value(s) s.

The current density j(s) = dI/ds is supplied as a callable. The enclosed current is its integral I(s) = ∫₀ˢ j ds’, and VMEC normalises so that I(1) corresponds to curtor. This function returns the dimensionless, normalised shape

Î(s) = (∫₀ˢ j ds’) / (∫₀¹ j ds’)

which runs from 0 at the axis to 1 at the plasma boundary. Dividing by the total ∫₀¹ j is essential: the density is only a shape, and its overall scale is arbitrary (it is commonly — but not always — normalised so that ∫₀¹ j = 1). Omitting this division is what made the old volume-current path wrong by the factor ∫₀¹ j.

Args:
current_density: callable j(s) = dI/ds (e.g. from power_series

of the VMEC ac coefficients).

s: scalar or 1-D array of normalised toroidal-flux values.

Returns:

Î(s) with the same shape as s (scalar in → scalar out). Returns all zeros when ∫₀¹ j ds == 0 (e.g. a zero current profile).

spectre.converters.vmec_helpers.distribute_current_to_interfaces(current_density, curtor, tfluxes, inface_current=False)

Map a continuous VMEC current profile onto SPECTRE interfaces.

Warning

Provisional / not final. How to map a continuous (VMEC) current profile onto the discrete interfaces of an MRxMHD (SPEC/SPECTRE) equilibrium is a fundamental open question. The control-volume scheme below is one physically-motivated choice; expect it to change.

inface_current = False (default) -> all current carried in the volumes; there is no net toroidal current on the interfaces. inface_current = True -> all current carried on the interfaces (isurf).

  • Volume current (inface_current=False) ivolume[i] = (−curtor·μ₀) · Î(sᵢ) — the cumulative enclosed current at interface i (Î from normalized_enclosed_current()); isurf is then zero.

  • Sheet current (inface_current=True) isurf[i] = (−curtor·μ₀) · (Î(eᵢ) − Î(eᵢ₋₁)) — the current density integrated over the control volume owned by interface i, whose edges sit at the midpoints to its neighbours:

    e₀ = 0                     # magnetic axis
eᵢ = (sᵢ + sᵢ₊₁)/2         # midpoint to the next interface
e_N = s_N                  # outer edge of the last cell = last interface

    So each interface collects current from halfway to its inner neighbour to halfway to its outer neighbour — an average of inside and outside, not just the current of the volume just inside it. For equidistant interfaces this reduces to the spot-evaluation j(sᵢ)·Δ; for non-equidistant interfaces a sparsely-spaced interface automatically owns a wider band and carries more of the surrounding current.

Current is conserved: the control-volume increments telescope to Î(s_N), so ivolume[-1] + sum(isurf) = (−curtor·μ₀)·Î(s_N) in either mode. When the last interface is the true plasma boundary (s_N = 1) this is the full −curtor·μ₀; when tfluxes[-1] < 1 the total is reduced by the appropriate amount, because the current flowing beyond the (interior) SPECTRE boundary is not enclosed.

The last interface carries K_N = (−curtor·μ₀)(Î(s_N) − Î(e_{N-1})) — the current from the last midpoint out to the boundary. In a fixed-boundary run SPECTRE ignores Isurf(Mvol), but this value is preserved so that fix2free() carries it over: when the plasma boundary is promoted to an interior interface in the free-boundary equilibrium, the slot becomes physical.

Args:
current_density: callable j(s) = dI/ds (e.g. from power_series

of the VMEC ac coefficients).

curtor: total toroidal current from VMEC [A]. tfluxes: 1-D array of interface toroidal-flux values in (0, 1]. inface_current: if False (default), all current is distributed as

volume current (ivolume); if True, all current is placed on the interfaces as sheet current (isurf).

Returns:

Tuple (ivolume, isurf) of 1-D arrays of length len(tfluxes), both in μ₀·A.

spectre.converters.vmec_helpers.vmec_boundary_to_fourier_coeffs(boundary, angle_fix=True)

Extract RBC and ZBS Fourier coefficients from a simsopt boundary surface.

A simsopt SurfaceRZFourier stores the full (m, n) mode spectrum including negative-m modes. We iterate over the positive-m half (num_modes = len(boundary.m) // 2 + 1) and read each cosine (R) and sine (Z) coefficient via the surface’s accessor methods.

Args:

boundary: simsopt SurfaceRZFourier object (e.g. vmec.boundary). angle_fix: if True (default), apply the θ→-θ handedness correction

for m>0 modes: the toroidal index is negated (n→-n) and the ZBS coefficient is negated. This pre-applies what lchangeangle=True would do at Fortran runtime, so the output TOML needs no lchangeangle flag. Set to False to keep the raw VMEC convention.

Returns:

Tuple (rbc, zbs) where each is a dict keyed by '(m, n)' strings mapping to the corresponding Fourier coefficient.

spectre.converters.vmec_helpers.compute_interface_guesses(wout, tfluxes, mpol, ntor)

Interpolate VMEC wout surfaces at each interface position.

Reads rmnc/zmns from the simsopt wout object, interpolates cubically on the VMEC radial grid, and evaluates at each value in tfluxes. The θ→-θ angle fix is applied (n→-n for m>0, ZBS negated) so the returned coefficients are in SPECTRE’s right-handed coordinate convention — identical to what lchangeangle=True would produce at Fortran runtime.

Only handles stellarator-symmetric equilibria (lasym=False).

Args:

wout: simsopt wout object (vmec_obj.wout). tfluxes: 1-D array of interface toroidal-flux values in (0, 1]. mpol: maximum poloidal mode number to include in the output. ntor: maximum toroidal mode number (reduced, without nfp) to include.

Returns:

dict suitable for PhysicsParameters.allrzrz: {'interface_1': {'rbc': {...}, 'zbs': {...}, 'rbs': {}, 'zbc': {}}, ...}

class spectre.converters.vmec_helpers.ContinuousPhysics

SPECTRE-relevant physics extracted from a VMEC equilibrium.

This is the contract between a VMEC adapter and build_spectre_input(). Every adapter must populate all fields; the builder has no knowledge of where the data came from.

The profiles are stored as callables f(s) on normalised toroidal flux s ∈ [0, 1], not as backend-specific coefficient arrays. This is what keeps the builder backend-agnostic: each adapter decides how to evaluate a profile (power series, grid interpolation, spline, …) and hands over a plain function. Being a mutable dataclass, a caller can override a profile in place, e.g. physics.iota = lambda s: 0.9 + 0.1 * s.

Attributes:

curtor: total toroidal current [A]. phiedge: toroidal magnetic flux at the plasma boundary [Wb]. nfp: number of toroidal field periods. rbc: cosine boundary coefficients as {'(m, n)': value} dict. zbs: sine boundary coefficients as {'(m, n)': value} dict. profile_type: 'iota' for a rotational-transform constraint,

'current' for an integrated-current constraint.

pressure: callable pressure(s) returning the pressure in SPECTRE

units (mu0-scaled); see pressure_function().

iota: callable iota(s) returning the rotational transform; see

iota_function(). Only used when profile_type == 'iota'.

current_density: callable j(s) = dI/ds (current-profile shape; its

overall scale is irrelevant, it is normalised internally). Only used when profile_type == 'current'.

curtor: float
phiedge: float
nfp: int
rbc: dict
zbs: dict
profile_type: str
pressure: collections.abc.Callable
iota: collections.abc.Callable
current_density: collections.abc.Callable
spectre.converters.vmec_helpers.build_spectre_input(physics, tfluxes, mpol, ntor, inface_current=False, interface_guesses=None)

Assemble a SPECTRE InputParameters object from VMEC physics.

This function has no dependency on any VMEC interface library; it operates purely on the ContinuousPhysics contract (profile callables + geometry).

Args:

physics: ContinuousPhysics dataclass instance. tfluxes: 1-D array of interface toroidal-flux values in (0, 1]. mpol: poloidal mode resolution for SPECTRE. ntor: toroidal mode resolution for SPECTRE. inface_current: if False (default), the toroidal current is

distributed as volume current (ivolume); if True, it is placed entirely on the interfaces as sheet current (isurf). Only used when physics.profile_type == 'current'. See distribute_current_to_interfaces() for the (provisional) discretisation scheme.

interface_guesses: optional dict in allrzrz format returned by

compute_interface_guesses(). When provided, it is stored as the initial geometry guess and linitialize is set to 0 (use supplied surfaces). When None, linitialize=1 causes SPECTRE to linearly interpolate an initial guess from the boundary.

Returns:

Tuple (dict_pydantic, profiles) where profiles is a dict of intermediate arrays ('pressure', 'iota', 'ivolume', 'isurf') useful for plotting or inspection.

Raises:

RuntimeError: if profile_type is unrecognised.