vmecpp package

exception vmecpp.HotRestartMismatchError

Bases: ValueError

vmecpp.run(input, magnetic_field=None, *, max_threads=None, verbose=OutputMode.PROGRESS, restart_from=None, iteration_callback=None)

Run VMEC++ using the provided input. This is the main entrypoint for both fixed- and free-boundary calculations.

Parameters:

  • input (VmecInput) – a VmecInput instance, corresponding to the contents of a classic VMEC input file

  • magnetic_field (Union[MagneticFieldResponseTable, Callable[[ndarray, ndarray, ndarray], tuple[ndarray, ndarray, ndarray]], None]) – the external magnetic field of a free-boundary run, used instead of the mgrid file in the input: either a MagneticFieldResponseTable, which is weighted by input.extcur and interpolated like an mgrid, or a MagneticFieldFunction b_r, b_phi, b_z = f(r, phi, z) evaluated at the plasma boundary points.

  • max_threads (Optional[int]) – maximum number of threads that VMEC++ should spawn. The actual number might still be lower that this in case there are too few flux surfaces to keep these many threads busy. If None, a number of threads equal to the number of logical cores is used.

  • verbose (bool | int | OutputMode) – controls the output format. Accepts bool for backward compatibility: 0. silent, no logging (False) 1. legacy table output (True) 2. animated progress bar for TTY enabled terminals 3. animated progress bar for non-TTY outputs (e.g. Jupyter)

  • restart_from (Union[VmecOutput, VmecWOut, None]) – if present, VMEC++ is initialized using the converged equilibrium from the provided VmecOutput or VmecWOut (e.g. loaded with VmecWOut.from_wout_file). This can dramatically decrease the number of iterations to convergence when running VMEC++ on a configuration that is very similar to the restart_from equilibrium. If input.mpol/input.ntor is a sequence (see below), this is used to hot-restart only the first continuation step; later steps always hot-restart from the previous one.

  • iteration_callback (Optional[Callable[[SolverState], bool | None]]) – called once per force iteration with a SolverState of the state just reached, after every thread has finished the step. Returning False stops the run, which then returns the outputs of that state with wout.ier_flag reporting no convergence; returning None or True continues. An exception raised inside the callback stops the run and propagates, and a return value that does not convert to a bool stops the run with a TypeError. With ncurr = 1 the callback may change the values of state.curr_h, the enclosed current the following force evaluations prescribe.

Return type:

VmecOutput

If input.mpol and/or input.ntor is a sequence rather than a plain int, run performs continuation in Fourier resolution: each entry pairs with the corresponding input.ns_array entry (a scalar mpol/ntor broadcasts to every step), and each step is solved in turn, hot-restarting from the previous step’s solution interpolated to the new resolution (see interpolate_solution).

Example

>>> import vmecpp
>>> path = "examples/data/solovev.json"
>>> vmec_input = vmecpp.VmecInput.from_file(path)
>>> output = vmecpp.run(vmec_input, verbose=False, max_threads=1)
>>> round(output.wout.b0, 10) # Exact value may differ by C library
0.2033313711
vmecpp.interpolate_solution(source, target_input)

Interpolate a converged solution onto the resolution of target_input.

Parameters:

  • source (VmecOutput) – a converged VmecOutput.

  • target_input (VmecInput) – the input for the next continuation step. Its mpol, ntor and (single-element) ns_array define the target resolution, and its nfp must match the source.

Return type:

VmecOutput

Returns:

A VmecOutput whose wout holds the source geometry interpolated to the target resolution, suitable as restart_from for vmecpp.run() with target_input.

vmecpp.rescale(output, b_scale, r_scale, scale_pressure=True, magnetic_field=None)
Parameters:
Return type:

VmecOutput

class vmecpp.VmecInput(**data)

Bases: BaseModelWithNumpy

The input to a VMEC++ run. Contains settings as well as the definition of the plasma boundary.

Python equivalent of a VMEC++ JSON input file or a classic INDATA file (e.g. “input.best”).

Deserialize from JSON and serialize to JSON using the usual pydantic methods: model_validate_json and model_dump_json.

Parameters:

data (Any)

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'ser_json_inf_nan': 'strings'}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

lasym: bool

Flag to indicate non-stellarator-symmetry.

  • False, assumes stellarator symmetry (only cosine/sine coefficients used).

  • True, (currently unsupported) allows for non-stellarator-symmetric terms.

nfp: int

Number of toroidal field periods (=1 for Tokamak)

mpol: MpolNtorField

Number of poloidal Fourier harmonics; m = 0, 1, …, (mpol-1).

May also be a sequence of ints, with one entry per ns_array step (a scalar broadcasts to every step), to request continuation in Fourier resolution: vmecpp.run() then solves each step in turn, hot-restarting from the previous step’s solution interpolated to the new resolution (see interpolate_solution()). The boundary coefficients (rbc, zbs, …) are always defined at the final (largest-index) entry’s resolution.

ntor: MpolNtorField

Number of toroidal Fourier harmonics; n = -ntor, -ntor+1, …, -1, 0, 1, …, ntor-1, ntor.

May be a sequence of ints, analogous to mpol; see its docstring.

keys()

The input field names, as for the dict-like INDATA objects of SIMSOPT.

Return type:

list[str]

property mpol_max: int

The final mpol resolution, if a multigrid sequence is used.

property ntor_max: int

The final ntor resolution, if a multigrid sequence is used.

mpol_geometry: int

Optional reduced poloidal resolution for the geometry (R, Z).

If in [1, mpol), R/Z modes with m >= mpol_geometry are held fixed while lambda keeps the full mpol. < 0 (default) means geometry uses mpol.

ntor_geometry: int

Optional reduced toroidal resolution for the geometry (R, Z).

If in [0, ntor), R/Z modes with n > ntor_geometry are held fixed while lambda keeps the full ntor. < 0 (default) means geometry uses ntor.

vacuum_mpol: int

Poloidal Fourier cutoff of the vacuum potential in a free-boundary run.

NESTOR expands the vacuum potential to mpol like the plasma; a value above mpol raises the potential’s cutoff alone. 0 (default) means the potential uses mpol.

vacuum_ntor: int

Toroidal Fourier cutoff of the vacuum potential in a free-boundary run.

A value above ntor raises the potential’s cutoff alone; nzeta must then be at least 2 * vacuum_ntor + 4. 0 (default) means the potential uses ntor.

ntheta: int

Number of poloidal grid points (ntheta >= 0).

Controls the poloidal resolution in real space. If 0, chosen automatically as minimally allowed. Must be at least 2*mpol + 6.

nzeta: int

Number of toroidal grid points (nzeta >= 0).

Controls the toroidal resolution in real space. If 0, chosen automatically as minimally allowed. Must be at least 2*ntor + 4. We typically use use phi as the convention for the toroidal angle, the name nzeta is due to beckwards compatibility.

ns_array: jt.Int[np.ndarray, 'num_grids']

Number of flux surfaces per multigrid step.

Each entry >= 3 and >= previous entry.

ftol_array: jt.Float[np.ndarray, 'num_grids']

Requested force tolerance for convergence per multigrid step.

niter_array: jt.Int[np.ndarray, 'num_grids']

Maximum number of iterations per multigrid step.

phiedge: float

Total enclosed toroidal magnetic flux in Vs == Wb.

  • In fixed-boundary, this determines the magnetic field strength.

  • In free-boundary, the magnetic field strength is given externally, so this determines cross-section area and volume of the plasma.

ncurr: typing.Literal[0, 1]

Select constraint on iota or enclosed toroidal current profiles.

  • 0: constrained-iota (rotational transform profile specified)

  • 1: constrained-current (toroidal current profile specified)

pmass_type: ProfileType

Parametrization of mass/pressure profile.

am: jt.Float[np.ndarray, 'am_len']

Mass/pressure profile coefficients.

Units: Pascals for pressure.

am_aux_s: jt.Float[np.ndarray, 'am_aux_len']

knot locations in s

Type:

Spline mass/pressure profile

am_aux_f: jt.Float[np.ndarray, 'am_aux_len']

values at knots

Type:

Spline mass/pressure profile

pres_scale: float

Global scaling factor for mass/pressure profile.

gamma: float

Adiabatic index \(\gamma\) (ratio of specific heats).

Specifying 0 implies that the pressure profile is specified. For all other values, the mass profile is specified.

spres_ped: float

Location of pressure pedestal in s.

Outside this radial location, pressure is constant.

piota_type: ProfileType

Parametrization of iota (rotational transform) profile.

ai: jt.Float[np.ndarray, 'ai_len']

Iota profile coefficients.

ai_aux_s: jt.Float[np.ndarray, 'ai_aux_len']

knot locations in s

Type:

Spline iota profile

ai_aux_f: jt.Float[np.ndarray, 'ai_aux_len']

values at knots

Type:

Spline iota profile

pcurr_type: ProfileType

Parametrization of toroidal current profile.

ac: jt.Float[np.ndarray, 'ac_len']

Enclosed toroidal current profile coefficients.

ac_aux_s: jt.Float[np.ndarray, 'ac_aux_len']

knot locations in s

Type:

Spline toroidal current profile

ac_aux_f: jt.Float[np.ndarray, 'ac_aux_len']

values at knots

Type:

Spline toroidal current profile

curtor: float

Net toroidal current in A.

The toroidal current profile is scaled to yield this total.

bloat: float

Bloating factor (for constrained toroidal current)

lfreeb: bool

Flag to indicate free-boundary.

If True, run in free-boundary mode; if False, fixed-boundary.

mgrid_file: typing.Annotated[str, pydantic.Field(max_length=200)]

Full path for vacuum Green’s function data.

NetCDF MGRID file with magnetic field response factors for external coils.

extcur: jt.Float[np.ndarray, 'ext_current']

Coil currents in A.

nvacskip: int

Number of iterations between full vacuum calculations.

signgs: int

-1 for the left-handed system of Fortran VMEC, +1 for a right-handed one.

The input boundary is flipped in theta to match it.

Type:

Sign of the Jacobian of the (s, theta, zeta) coordinates

free_boundary_method: typing.Annotated[FreeBoundaryMethod, pydantic.BeforeValidator(_validate_free_boundary_method), pydantic.Field()]

Method for handling free-boundary conditions.

iteration_style: typing.Annotated[IterationStyle, pydantic.BeforeValidator(_validate_iteration_style), pydantic.Field()]

Time-step / restart control scheme for the equilibrium iteration ("vmec_8_52" or "parvmec").

nstep: int

Printout interval at which convergence progress is logged.

aphi: jt.Float[np.ndarray, 'aphi_len']

Radial flux zoning profile coefficients.

delt: float

Initial value for artificial time step in iterative solver.

tcon0: float

Constraint force scaling factor for ns –> 0.

lgiveup: bool

Abandon the whole multigrid sequence when a step ends with any residual still above fgiveup times its tolerance, rather than carrying a state that far out onto a finer grid.

fgiveup: float

Multiple of ftol_array a step’s residuals must be under for the sequence to continue when lgiveup is set.

lforbal: bool

directly compute innermost flux surface geometry from radial force balance

Type:

Hack

lambda_preconditioner_scale: float

Scale of the lambda preconditioner, which multiplies the inverse of the diagonal lambda stiffness to turn the lambda force into the lambda step.

1.0 applies the undamped inverse, values below 1.0 damp the lambda step and values above 1.0 accelerate it. The default 0.5 is the damping of VMEC 8.52.

lambda_precondition_checkerboard_terms: float

Stiffness floor with which the lambda preconditioner treats the radial checkerboard (odd-even) mode of the m = 0 lambda, as a fraction of the stiffness of a smooth mode.

lambda enters B through its average over neighbouring surfaces, which cancels a (-1)^j checkerboard of the full-grid m = 0 lambda, so the energy and the force hardly see it; a diagonal lambda preconditioner converges it last, and it shows up as a sign-alternating <J.B> and DMerc. Non-negative values precondition the m = 0 lambda with its tridiagonal radial block plus this fraction of the diagonal, which enlarges the checkerboard step at most 1 / value times: 0 inverts the exact block. Negative values keep the diagonal preconditioner of VMEC 8.52, which the reference tests use. The converged equilibrium does not depend on it.

lbsubs: bool

If true, recompute the full-grid covariant B_s by solving radial force balance with regularized radial derivatives.

Helps with getting more physical currents for post-processing involving B_s like diagnostic reconstruction, DMerc, or toroidal current densities. https://proximafusion.github.io/vmecpp/tutorials/fixed_boundary.html

backup_evaluated_state: bool

If true, restart backups hold the state of the last force evaluation; if false, they hold the advanced state, as in educational_VMEC and the reference files.

always_fix_m1_gauge: bool

If true, the m = 1 gauge of R and Z is set from the boundary at the start of every multigrid step and held there; if false, it is held only in the first iteration of a multigrid step and while fsqz < 1e-6, as in VMEC 8.52.

return_outputs_even_if_not_converged: bool

If true, return a wout even if VMEC++ did not converge, instead of raising a RuntimeError.

This is intended for debugging purposes (e.g. inspecting how far the geometry got, or where the force residuals blew up) since the returned quantities are computed from whatever internal state the solver was in when it gave up, and can be arbitrarily unphysical. Always check wout.ier_flag and the residuals wout.fsqr, wout.fsqz, wout.fsql against wout.ftolv before interpreting any physical quantity in the output.

raxis_c: jt.Float[np.ndarray, 'ntor_plus_1']

Magnetic axis coefficients for R ~ cos(n*v); stellarator-symmetric.

At least 1 value required, up to n=ntor considered.

zaxis_s: jt.Float[np.ndarray, 'ntor_plus_1']

Magnetic axis coefficients for Z ~ sin(n*v); stellarator-symmetric.

Up to n=ntor considered; first entry (n=0) is ignored.

raxis_s: jt.Float[np.ndarray, 'ntor_plus_1'] | None

Magnetic axis coefficients for R ~ sin(n*v); non-stellarator-symmetric.

Up to n=ntor considered; first entry (n=0) is ignored. Only used if lasym=True.

zaxis_c: jt.Float[np.ndarray, 'ntor_plus_1'] | None

Magnetic axis coefficients for Z ~ cos(n*v); non-stellarator-symmetric.

Only used if lasym=True.

property raxis_cc
property raxis_cs
property zaxis_cc
property zaxis_cs
rbc: SerializableSparseCoefficientArray[jt.Float[NpOrJax, 'mpol two_ntor_plus_one']]

Boundary coefficients for R ~ cos(m*u - n*v); stellarator-symmetric

zbs: SerializableSparseCoefficientArray[jt.Float[NpOrJax, 'mpol two_ntor_plus_one']]

Boundary coefficients for Z ~ sin(m*u - n*v); stellarator-symmetric

rbs: SerializableSparseCoefficientArray[jt.Float[np.ndarray, 'mpol two_ntor_plus_one']] | None

Boundary coefficients for R ~ sin(m*u - n*v); non-stellarator-symmetric.

Only used if lasym=True.

zbc: SerializableSparseCoefficientArray[jt.Float[np.ndarray, 'mpol two_ntor_plus_one']] | None

Boundary coefficients for Z ~ cos(m*u - n*v); non-stellarator-symmetric.

Only used if lasym=True.

static resize_1d_axis_coeff(coeff, ntor_new)

Resizes a 1D magnetic-axis Fourier coefficient array to length ntor_new+1.

Arrays shorter than ntor_new+1 are zero-padded (the omitted trailing coefficients are implicitly zero). Arrays longer than ntor_new+1 are truncated.

Parameters:

  • coeff (Float[ndarray, 'ntor_plus_1']) – A 1D NumPy array of axis coefficients (length ntor+1).

  • ntor_new (int) – The new number of toroidal modes.

Return type:

Float[ndarray, 'ntor_new_plus_1']

Examples

>>> VmecInput.resize_1d_axis_coeff(np.array([1.0, 2.0]), ntor_new=3)
array([1., 2., 0., 0.])
>>> VmecInput.resize_1d_axis_coeff(np.array([1.0, 2.0, 3.0]), ntor_new=1)
array([1., 2.])
static resize_2d_coeff(coeff, mpol_new, ntor_new)

Resizes a 2D NumPy array representing Fourier coefficients, padding with zeros or truncating as needed.

Parameters:

  • coeff (Float[ndarray, 'mpol two_ntor_plus_one']) – A NumPy array of shape (mpol, 2 * ntor + 1).

  • mpol_new (int) – The new number of poloidal modes.

  • ntor_new (int) – The new number of toroidal modes.

Return type:

Float[ndarray, 'mpol_new two_ntor_new_plus_one']

Examples

>>> coeff = np.array([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10]])
>>> VmecInput.resize_2d_coeff(coeff, 3, 3)
array([[ 0.,  1.,  2.,  3.,  4.,  5.,  0.],
       [ 0.,  6.,  7.,  8.,  9., 10.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]])

>>> VmecInput.resize_2d_coeff(coeff, 1, 1)
array([[2., 3., 4.]])

>>> VmecInput.resize_2d_coeff(coeff, 4, 1)
array([[2., 3., 4.],
       [7., 8., 9.],
       [0., 0., 0.],
       [0., 0., 0.]])
resize(mpol_new, ntor_new)

Resample this input in place to a new (mpol, ntor) Fourier resolution.

An mpol/ntor continuation schedule keeps its length: the last entry becomes the new size and earlier entries are clamped to it.

Boundary coefficients are zero-padded or truncated to match, discarding higher modes with a warning; see resize_2d_coeff(). Axis coefficients are zero-padded or, when shrinking ntor, truncated with a warning.

Parameters:

  • mpol_new (int)

  • ntor_new (int)

Return type:

None

static from_file(input_file)

Build a VmecInput from either a VMEC++ JSON input file or a classic INDATA file.

Parameters:

input_file (str | Path)

Return type:

VmecInput

to_json(**kwargs)

Serialize the object to JSON.

Keyword Arguments:

**kwargs – Additional keyword arguments to forward to the model_dump_json method.

Return type:

str

save(output_path, **kwargs)
Parameters:

output_path (str | Path)

Return type:

None

class vmecpp.VmecOutput(**data)

Bases: BaseModelWithNumpy

Container for the full output of a VMEC run.

Parameters:

data (Any)

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'ser_json_inf_nan': 'strings'}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

input: VmecInput

The input to the VMEC run that produced this output.

jxbout: JxBOut

Python equivalent of VMEC’s “jxbout” file.

mercier: Mercier

Python equivalent of VMEC’s “mercier” file.

Contains radial profiles and stability criteria relevant for Mercier stability analysis, including the Mercier criterion and its decomposition into shear, well, current, and geodesic contributions. Also includes profiles of rotational transform, toroidal flux, pressure, and their derivatives.

threed1_volumetrics: Threed1Volumetrics

Python equivalent of VMEC’s volumetrics section in the “threed1” file.

Contains global and flux-surface-averaged quantities such as total and average pressure, poloidal and toroidal magnetic field energies, kinetic energy, and related integrals. Useful for postprocessing and global equilibrium characterization.

threed1_first_table: Threed1FirstTable

radial profiles of s, iota, currents, pressure, and force balance.

Type:

Python equivalent of the first table in VMEC’s “threed1” file

threed1_geometric_magnetic: Threed1GeometricAndMagneticQuantities

Python equivalent of the geometric and magnetic quantities in VMEC’s “threed1” file: global geometry, beta values, currents, and geometric profiles.

threed1_axis: Threed1AxisGeometry

Python equivalent of the magnetic-axis geometry in VMEC’s “threed1” file.

threed1_betas: Threed1Betas

Python equivalent of the beta values in VMEC’s “threed1” file.

threed1_shafranov_integrals: Threed1ShafranovIntegrals

Python equivalent of the Shafranov surface integrals in VMEC’s “threed1” file.

wout: VmecWOut

Python equivalent of VMEC’s “wout” file.

class vmecpp.VmecWOut(**data)

Bases: BaseModelWithNumpy

Python equivalent of a VMEC “wout file”.

VmecWOut exposes the layout that SIMSOPT expects. The save method produces a NetCDF file compatible with SIMSOPT/Fortran VMEC wout.nc.

Parameters:

data (Any)

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'allow', 'ser_json_inf_nan': 'strings', 'serialize_by_alias': False, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

input_extension: typing.Annotated[str, pydantic.Field(max_length=100)]

File extension of the input file.

ier_flag: int

Status code indicating success or problems during the VMEC++ run.

See the reason property for a human-readable description.

property reason: str
nfp: int

Number of toroidal field periods.

ns: int

Number of radial grid points (=number of flux surfaces).

mpol: int

Number of poloidal Fourier modes.

ntor: int

Number of toroidal Fourier modes.

mnmax: int

Number of Fourier coefficients for the state vector.

mnmax_nyq: int

Number of Fourier coefficients for the Nyquist-quantities.

lasym: typing.Annotated[bool, pydantic.PlainSerializer(int), pydantic.BeforeValidator(bool), pydantic.Field(alias='lasym__logical__')]

Flag indicating non-stellarator-symmetry.

Non-stellarator symmetric fields are only populated if this is True.

lfreeb: typing.Annotated[bool, pydantic.PlainSerializer(int), pydantic.BeforeValidator(bool), pydantic.Field(alias='lfreeb__logical__')]

Flag indicating free-boundary computation.

lrfp: typing.Annotated[bool, pydantic.PlainSerializer(int), pydantic.BeforeValidator(bool), pydantic.Field(alias='lrfp__logical__', default=False)]

Flag indicating reversed-field pinch configuration.

wb: FloatOrJax

volume integral of |B|^2/(2 mu0).

Type:

Magnetic energy

wp: FloatOrJax

volume integral of p.

Type:

Kinetic energy

rmax_surf: FloatOrJax

Maximum R on the plasma boundary over all grid points.

rmin_surf: FloatOrJax

Minimum R on the plasma boundary over all grid points.

zmax_surf: FloatOrJax

Maximum Z on the plasma boundary over all grid points.

aspect: FloatOrJax

Aspect ratio (major radius over minor radius) of the plasma boundary.

betapol: FloatOrJax

Poloidal plasma beta.

The ratio of the total thermal energy of the plasma to the total poloidal magnetic energy. \(\beta = W_{th} / W_{B_\theta} = \int p\, dV / \left( \int B_\theta^2 / (2 \mu_0)\, dV \right )\)

betator: FloatOrJax

Toroidal plasma beta.

The ratio of the total thermal energy of the plasma to the total toroidal magnetic energy. \(\beta = W_{th} / W_{B_\phi} = \int p\, dV / \left( \int B_\phi^2 / (2 \mu_0)\, dV \right )\)

betaxis: FloatOrJax

Plasma beta on the magnetic axis.

b0: FloatOrJax

Toroidal magnetic flux density from poloidal current and magnetic axis position at phi=0.

rbtor0: FloatOrJax

Poloidal ribbon current at the axis.

rbtor: FloatOrJax

Poloidal ribbon current at the plasma boundary.

IonLarmor: FloatOrJax

Larmor radius of plasma ions.

ctor: FloatOrJax

Net toroidal plasma current.

Aminor_p: FloatOrJax

Minor radius of the plasma.

Rmajor_p: FloatOrJax

Major radius of the plasma.

volume: typing.Annotated[FloatOrJax, pydantic.Field(alias='volume_p')]

Plasma volume.

fsqr: FloatOrJax

Invariant force residual of the force on R at end of the run.

fsqz: FloatOrJax

Invariant force residual of the force on Z at end of the run.

fsql: FloatOrJax

Invariant force residual of the force on lambda at end of the run.

ftolv: FloatOrJax

Force tolerance value used to determine convergence.

itfsq: int

Number of force-balance iterations after which the run terminated.

phipf: jt.Float[NpOrJax, 'n_surfaces']

Radial derivative of enclosed toroidal magnetic flux phi' on the full- grid.

chipf: jt.Float[NpOrJax, 'n_surfaces']

Radial derivative of enclosed poloidal magnetic flux chi' on the full- grid.

jcuru: jt.Float[NpOrJax, 'n_surfaces']

Radial derivative of enclosed poloidal current on full-grid.

jcurv: jt.Float[NpOrJax, 'n_surfaces']

Radial derivative of enclosed toroidal current on full-grid.

fsqt: jt.Float[NpOrJax, 'time']

Evolution of the total force residual along the run.

This is the sum of force_residual_r, force_residual_z, and force_residual_lambda.

force_residual_r: jt.Float[NpOrJax, 'time']

Evolution of the r radial force residual along the run.

force_residual_z: jt.Float[NpOrJax, 'time']

Evolution of the z vertical force residual along the run.

force_residual_lambda: jt.Float[NpOrJax, 'time']

Evolution of the lambda force residual along the run.

delbsq: jt.Float[NpOrJax, 'time']

Evolution of the force residual at the vacuum boundary along the run.

lcfs_force: jt.Float[NpOrJax, 'time']

Evolution of the normalized squared R and Z force residuals on the LCFS.

restart_reason_timetrace: typing.Annotated[jt.Int[NpOrJax, 'time'], pydantic.Field(alias='restart_reasons'), pydantic.BeforeValidator(lambda x: np.array(x).astype(np.int64))]

Internal restart reasons at each step along the run. (debugging quantity).

Use the restart_reasons field to access a more readable enum version of this instead of integer status codes.

wdot: jt.Float[NpOrJax, 'time']

Evolution of the MHD energy decay along the run.

jdotb: jt.Float[NpOrJax, 'n_surfaces']

Flux-surface-averaged \(\langle j \cdot B \rangle\) on full-grid.

bdotb: jt.Float[NpOrJax, 'n_surfaces']

Flux-surface-averaged \(\langle B \cdot B \rangle\) on full-grid.

bdotgradv: jt.Float[NpOrJax, 'n_surfaces']

Flux-surface-averaged toroidal magnetic field component \(B \cdot \nabla v\) on full-grid.

DMerc: jt.Float[NpOrJax, 'n_surfaces']

Full Mercier stability criterion on the full-grid.

equif: jt.Float[NpOrJax, 'n_surfaces']

Radial force balance residual on full-grid.

potvac: jt.Float[NpOrJax, 'mnpd'] | None

Fourier coefficients of Nestor’s scalar magnetic potential, sin(mu - nv) followed by cos(mu - nv). None for a fixed-boundary run.

xmpot: SerializeIntAsFloat[jt.Int[NpOrJax, 'mn_mode_pot']] | None

Poloidal mode numbers m for the potvac coefficients.

None for a fixed-boundary run.

xnpot: SerializeIntAsFloat[jt.Int[NpOrJax, 'mn_mode_pot']] | None

Toroidal mode numbers times number of toroidal field periods n * nfp for the potvac coefficients.

None for a fixed-boundary run.

xm: SerializeIntAsFloat[jt.Int[NpOrJax, 'mn_mode']]

Poloidal mode numbers m for the Fourier coefficients in the state vector.

xn: SerializeIntAsFloat[jt.Int[NpOrJax, 'mn_mode']]

Toroidal mode numbers times number of toroidal field periods n * nfp for the Fourier coefficients in the state vector.

xm_nyq: SerializeIntAsFloat[jt.Int[NpOrJax, 'mn_mode_nyq']]

Poloidal mode numbers m for the Fourier coefficients in the Nyquist- quantities.

xn_nyq: SerializeIntAsFloat[jt.Int[NpOrJax, 'mn_mode_nyq']]

Toroidal mode numbers times number of toroidal field periods n * nfp for the Fourier coefficients in the Nyquist-quantities.

mass: jt.Float[NpOrJax, 'n_surfaces']

Plasma mass profile m on half-grid.

buco: jt.Float[NpOrJax, 'n_surfaces']

Profile of enclosed toroidal current I on half-grid.

bvco: jt.Float[NpOrJax, 'n_surfaces']

Profile of enclosed poloidal ribbon current G on half-grid.

phips: jt.Float[NpOrJax, 'n_surfaces']

Radial derivative of enclosed toroidal magnetic flux phi' on the half- grid.

bmnc: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces']

Fourier coefficients (cos) of the magnetic field strength |B| on the half- grid.

gmnc: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces']

Fourier coefficients (cos) of the Jacobian \(\sqrt{g}\) on the half-grid.

bsubumnc: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces']

Fourier coefficients (cos) of the covariant magnetic field component \(B_{\theta}\) on the half-grid.

bsubvmnc: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces']

Fourier coefficients (cos) of the covariant magnetic field component \(B_{\phi}\) on the half-grid.

bsubsmns: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces']

Fourier coefficients (sin) of the covariant magnetic field component \(B_{s}\) on the full-grid; Fortran VMEC 8.52 writes the half-grid array.

bsupumnc: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces']

Fourier coefficients (cos) of the contravariant magnetic field component \(B^{\theta}\) on the half-grid.

bsupvmnc: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces']

Fourier coefficients (cos) of the contravariant magnetic field component \(B^{\phi}\) on the half-grid.

currumnc: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces']

Fourier coefficients (cos) of \(\sqrt{g} J^{\theta}\) on the full-grid.

Evaluate the cosine series with xm_nyq * theta - xn_nyq * zeta, add the sine series from currumns when present, then divide by the Jacobian \(\sqrt{g}\) at the same points to obtain \(J^{\theta}\). xn_nyq already includes nfp. These coefficients include the factor \(1/\mu_0\); do not divide by \(\mu_0\) again. Prefer these exported currents when computing current-dependent quantities.

At an interior full-grid index j, the matching Jacobian is obtained by averaging half-grid Fourier coefficients at columns j and j + 1 of gmnc and, when present, gmns. Axis and boundary currents are extrapolated values.

Consumers rebuilding currents from covariant B must use the same radial derivative as VMEC++. With lbsubs=True, the full-grid B_s is paired with the sqrt(s)-regularized odd-m derivative of B_u and B_v. For either component Q, write Q = sqrt(s) * q for odd m and evaluate at interior surface j:

h = 1 / (ns - 1)
s = j * h
q_inner = Q[j] / sqrt((j - 0.5) * h)
q_outer = Q[j + 1] / sqrt((j + 0.5) * h)
dQ_ds = sqrt(s) * (q_outer - q_inner) / h
dQ_ds += (q_inner + q_outer) / (4 * sqrt(s))

Here Q is one Fourier coefficient across the stored half-grid columns; apply this to both sine and cosine coefficients. Even m uses the ordinary difference (Q[j + 1] - Q[j]) / h. Angular derivatives of B_s use the Fourier basis. The same rule applies to other calculations combining these radial derivatives with B_s. Plain odd-m differences break the cancellation at the axis and can produce a spurious current growing as \(h^{-1/2}\) under radial refinement.

currvmnc: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces']

Fourier coefficients (cos) of \(\sqrt{g} J^{\zeta}\) on the full-grid.

Reconstruct with xm_nyq, xn_nyq, and the sine partner currvmns, then divide by the full-grid Jacobian. The physical toroidal current density is \(J_{\phi} = R J^{\zeta}\) in A/m^2. The factor \(1/\mu_0\) is already included. See currumnc for Jacobian interpolation and the required sqrt(s)-regularized odd-m derivative when rebuilding currents from B_s, B_u, and B_v.

rmnc: jt.Float[NpOrJax, 'mn_mode n_surfaces']

Fourier coefficients (cos) for R of the geometry of the flux surfaces on the full- grid.

zmns: jt.Float[NpOrJax, 'mn_mode n_surfaces']

Fourier coefficients (sin) for Z of the geometry of the flux surfaces on the full- grid.

lmns: jt.Float[NpOrJax, 'mn_mode n_surfaces']

Fourier coefficients (sin) for lambda stream function on the half-grid.

lmns_full: jt.Float[NpOrJax, 'mn_mode n_surfaces']

Fourier coefficients (sin) for lambda stream function on the full-grid.

This quantity is VMEC++ specific and required for hot-restart to work properly. We store it with the Fortran convention for the order of the dimensions for consistency with lmns.

rmns: jt.Float[NpOrJax, 'mn_mode n_surfaces'] | None

Fourier coefficients (sin) for R of the geometry of the flux surfaces on the full-grid; non-stellarator-symmetric.

zmnc: jt.Float[NpOrJax, 'mn_mode n_surfaces'] | None

Fourier coefficients (cos) for Z of the geometry of the flux surfaces on the full-grid; non-stellarator-symmetric.

lmnc: jt.Float[NpOrJax, 'mn_mode n_surfaces'] | None

Fourier coefficients (cos) for lambda stream function on the half-grid; non- stellarator-symmetric.

lmnc_full: jt.Float[NpOrJax, 'mn_mode n_surfaces'] | None

Fourier coefficients (cos) for lambda stream function on the full-grid; non- stellarator-symmetric.

This quantity is VMEC++ specific and required for hot-restart to work properly. We store it with the Fortran convention for the order of the dimensions for consistency with lmnc. Only populated when lasym=True (non-stellarator-symmetric configurations).

gmns: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces'] | None

Fourier coefficients (sin) of the Jacobian \(\sqrt{g}\) on the half-grid; non-stellarator-symmetric.

bmns: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces'] | None

Fourier coefficients (sin) of the magnetic field strength |B| on the half- grid; non-stellarator-symmetric.

bsubumns: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces'] | None

Fourier coefficients (sin) of the covariant magnetic field component \(B_{\theta}\) on the half-grid; non-stellarator-symmetric.

bsubvmns: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces'] | None

Fourier coefficients (sin) of the covariant magnetic field component \(B_{\phi}\) on the half-grid; non-stellarator-symmetric.

bsubsmnc: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces'] | None

Fourier coefficients (cos) of the covariant magnetic field component \(B_{s}\) on the full-grid; non-stellarator-symmetric.

bsupumns: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces'] | None

Fourier coefficients (sin) of the contravariant magnetic field component \(B^{\theta}\) on the half-grid; non-stellarator-symmetric.

bsupvmns: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces'] | None

Fourier coefficients (sin) of the contravariant magnetic field component \(B^{\phi}\) on the half-grid; non-stellarator-symmetric.

currumns: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces'] | None

Fourier coefficients (sin) of \(\sqrt{g} J^{\theta}\) on the full-grid; non- stellarator-symmetric.

Use together with currumnc as described there.

currvmns: jt.Float[NpOrJax, 'mn_mode_nyq n_surfaces'] | None

Fourier coefficients (sin) of \(\sqrt{g} J^{\zeta}\) on the full-grid; non- stellarator-symmetric.

Use together with currvmnc as described there.

pcurr_type: ProfileType

Parametrization of toroidal current profile (copied from input).

pmass_type: ProfileType

Parametrization of mass/pressure profile (copied from input).

piota_type: ProfileType

Parametrization of iota profile (copied from input).

am: ProfileCoeffType[jt.Float[NpOrJax, '_preset']]

Mass/pressure profile coefficients (copied from input).

ac: ProfileCoeffType[jt.Float[NpOrJax, '_preset']]

Enclosed toroidal current profile coefficients (copied from input).

ai: ProfileCoeffType[jt.Float[NpOrJax, '_preset']]

Iota profile coefficients (copied from input).

am_aux_s: AuxSType[jt.Float[NpOrJax, '_ndfmax']]

knot locations in s (copied from input).

Type:

Spline mass/pressure profile

am_aux_f: AuxFType[jt.Float[NpOrJax, '_ndfmax']]

values at knots (copied from input).

Type:

Spline mass/pressure profile

ac_aux_s: AuxSType[jt.Float[NpOrJax, '_ndfmax']]

knot locations in s (copied from input).

Type:

Spline toroidal current profile

ac_aux_f: AuxFType[jt.Float[NpOrJax, '_ndfmax']]

values at knots (copied from input).

Type:

Spline toroidal current profile

ai_aux_s: AuxSType[jt.Float[NpOrJax, '_ndfmax']]

knot locations in s (copied from input).

Type:

Spline iota profile

ai_aux_f: AuxFType[jt.Float[NpOrJax, '_ndfmax']]

values at knots (copied from input).

Type:

Spline iota profile

gamma: FloatOrJax

Adiabatic index \(\gamma\) (copied from input).

mgrid_file: typing.Annotated[str, pydantic.Field(max_length=200)]

Full path for vacuum Green’s function data (copied from input).

nextcur: int

Number of external coil currents.

extcur: typing.Annotated[jt.Float[NpOrJax, 'ext_current'], pydantic.BeforeValidator(lambda x: x if np.shape(x) != () else np.array([])), pydantic.WrapSerializer(lambda x, handler, _: netCDF4.default_fillvals['f8'] if isinstance(x, (np.ndarray, list)) and np.shape(x) == (0,) else handler(x))]

Coil currents in A.

for free-boundary runs, extcur has shape (nextcur,) for fixed-boundary it is a scalar float extcur=nan

mgrid_mode: MgridModeType

Indicates if the mgrid file was normalized to unit currents (“S”) or not (“R”).

iotas: jt.Float[NpOrJax, 'n_surfaces']

Rotational transform \(\iota\) on the half-grid.

iotaf: jt.Float[NpOrJax, 'n_surfaces']

Rotational transform \(\iota\) on the full-grid.

betatotal: FloatOrJax

Total plasma beta.

The ratio of the total thermal energy of the plasma to the total magnetic energy.

\(\beta = W_{th} / W_B = \int p\, dV / \left( \int B^2 / (2 \mu_0)\, dV \right )\)

raxis_cc: jt.Float[NpOrJax, 'ntor_plus_1']

Fourier coefficients of \(R(phi)\) of the magnetic axis geometry.

zaxis_cs: jt.Float[NpOrJax, 'ntor_plus_1']

Fourier coefficients of \(Z(phi)\) of the magnetic axis geometry.

raxis_cs: jt.Float[NpOrJax, 'ntor_plus_1'] | None

Fourier coefficients of \(R(phi)\) of the magnetic axis geometry; non- stellarator-symmetric.

zaxis_cc: jt.Float[NpOrJax, 'ntor_plus_1'] | None

Fourier coefficients of \(Z(phi)\) of the magnetic axis geometry; non- stellarator-symmetric.

vp: jt.Float[NpOrJax, 'n_surfaces']

Differential volume \(V' = \frac{\partial V}{\partial s}\) on half-grid.

presf: jt.Float[NpOrJax, 'n_surfaces']

Kinetic pressure p on the full-grid.

pres: jt.Float[NpOrJax, 'n_surfaces']

Kinetic pressure p on the half-grid.

phi: jt.Float[NpOrJax, 'n_surfaces']

Enclosed toroidal magnetic flux \(\phi\) on the full-grid.

signgs: int

Sign of the Jacobian of the coordinate transform between flux coordinates and cylindrical coordinates.

volavgB: FloatOrJax

Volume-averaged magnetic field strength.

q_factor: jt.Float[NpOrJax, 'n_surfaces']

Safety factor \(q = 1/\iota\) on the full-grid.

chi: jt.Float[NpOrJax, 'n_surfaces']

Enclosed poloidal magnetic flux \(\chi\) on the full-grid.

specw: jt.Float[NpOrJax, 'n_surfaces']

Spectral width M on the full-grid.

over_r: jt.Float[NpOrJax, 'n_surfaces']

<\tau / R> / V' on half-grid.

\(\left\langle \frac{\tau}{R} \right\rangle / V'\)

DShear: jt.Float[NpOrJax, 'n_surfaces']

Mercier stability criterion contribution due to magnetic shear.

DWell: jt.Float[NpOrJax, 'n_surfaces']

Mercier stability criterion contribution due to magnetic well.

DCurr: jt.Float[NpOrJax, 'n_surfaces']

Mercier stability criterion contribution due to plasma currents.

DGeod: jt.Float[NpOrJax, 'n_surfaces']

Mercier stability criterion contribution due to geodesic curvature.

niter: int

Number of force-balance iterations taken to converge.

beta_vol: jt.Float[NpOrJax, 'n_surfaces']

Flux-surface averaged plasma beta on half-grid.

version_: FloatOrJax

Version number of VMEC, that this VMEC++ wout file is compatible with.

Some codes change how they interpret values in the wout file depending on this number. (E.g. COBRAVMEC checks if >6 or not)

property volume_p

The attribute is called volume_p in the Fortran wout file, while simsopt.mhd.Vmec.wout uses volume.

We expose both.

property lasym__logical__

This is how the attribute is called in the Fortran wout file.

property lfreeb__logical__

This is how the attribute is called in the Fortran wout file.

property lrfp__logical__

This is how the attribute is called in the Fortran wout file.

property restart_reasons: list[tuple[int, RestartReason]]

Get the restart reasons as a list of tuples.

Each tuple contains the iteration number and the reason for the restart.

save(out_path)

Save contents in NetCDF3 format, e.g. wout.nc.

This is the format used by Fortran VMEC implementations and the one expected by SIMSOPT.

Parameters:

out_path (str | Path)

Return type:

None

static from_wout_file(wout_filename)

Load wout contents in NetCDF format.

This is the format used by Fortran VMEC implementations and the one expected by SIMSOPT. We allow for additional attributes to be present in the file, for compatibility with wouf files from other VMEC versions, but require at least the fields produced by VMEC++.

Parameters:

wout_filename (str | Path)

Return type:

VmecWOut

class vmecpp.JxBOut(**data)

Bases: BaseModelWithNumpy

Parameters:

data (Any)

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'ser_json_inf_nan': 'strings'}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

itheta: jt.Float[np.ndarray, 'num_full nZnT']

Poloidal surface current.

\(itheta = (\frac{\partial B_s}{\partial \Phi} - \frac{\partial B_\phi}{\partial s}) / \mu_0\)

izeta: jt.Float[np.ndarray, 'num_full nZnT']

Toroidal surface current.

\(izeta = (-\frac{\partial B_s}{\partial \Theta} + \frac{\partial B_\theta}{\partial s}) / \mu_0\)

bdotk: jt.Float[np.ndarray, 'num_full nZnT']
amaxfor: jt.Float[np.ndarray, 'n_surfaces']

100 times the maximum value of the real space force residual on each radial surface.

aminfor: jt.Float[np.ndarray, 'n_surfaces']

100 times the minimum value of the real space force residual on each radial surface.

avforce: jt.Float[np.ndarray, 'n_surfaces']

Average force residual on each radial surface.

pprim: jt.Float[np.ndarray, 'n_surfaces']

Radial derivative of the pressure profile.

jdotb: jt.Float[np.ndarray, 'n_surfaces']

Flux-surface-averaged \(\langle j \cdot B \rangle\) on full-grid.

bdotb: jt.Float[np.ndarray, 'n_surfaces']

Flux-surface-averaged \(\langle B \cdot B \rangle\) on full-grid.

bdotgradv: jt.Float[np.ndarray, 'n_surfaces']
jpar2: jt.Float[np.ndarray, 'n_surfaces']

langle j_{||}^2 rangle.

Type:

Flux-surface-averaged squared parallel current density

Type:

math

jperp2: jt.Float[np.ndarray, 'n_surfaces']

langle j_{perp}^2 rangle.

Type:

Flux-surface-averaged squared perpendicular current density

Type:

math

phin: jt.Float[np.ndarray, 'n_surfaces']

Normalized, enclosed toroidal flux at each radial surface.

phin = toroidal_flux/toroidal_flux[-1]

jsupu3: jt.Float[np.ndarray, 'num_full nZnT']

Contravariant current density component j^u on the full grid.

\(j^u = itheta/V'\)

jsupv3: jt.Float[np.ndarray, 'num_full nZnT']

Contravariant current density component j^v on the full grid.

\(j^u = izeta/V'\)

jsups3: jt.Float[np.ndarray, 'num_half nZnT']

Contravariant current density component \(j^s\) on the half grid.

\(j^s = \frac{\partial B_\theta}{\partial \Phi} - \frac{\partial B_\phi}{\partial \Theta}{\mu_0 V'}\)

bsupu3: jt.Float[np.ndarray, 'num_full nZnT']
bsupv3: jt.Float[np.ndarray, 'num_full nZnT']
jcrossb: jt.Float[np.ndarray, 'num_full nZnT']

Magnitude of \(j \times B\) at each grid point.

jxb_gradp: jt.Float[np.ndarray, 'num_full nZnT']

Dot product of \(j \times B\) and \(\nabla p\) at each grid point.

jdotb_sqrtg: jt.Float[np.ndarray, 'num_full nZnT']

Product of \(j \cdot B\) and \(\sqrt{g}\) at each grid point.

sqrtg3: jt.Float[np.ndarray, 'num_full nZnT']

Jacobian determinant \(\sqrt{g}\) at each grid point.

bsubu3: jt.Float[np.ndarray, 'num_half nZnT']
bsubv3: jt.Float[np.ndarray, 'num_half nZnT']
bsubs3: jt.Float[np.ndarray, 'num_full nZnT']
class vmecpp.Mercier(**data)

Bases: BaseModelWithNumpy

Parameters:

data (Any)

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'ser_json_inf_nan': 'strings'}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

s: jt.Float[np.ndarray, 'n_surfaces']

Normalized toroidal flux coordinate s.

toroidal_flux: jt.Float[np.ndarray, 'n_surfaces']

Enclosed toroidal magnetic flux phi.

iota: jt.Float[np.ndarray, 'n_surfaces']

Rotational transform iota.

shear: jt.Float[np.ndarray, 'n_surfaces']

Magnetic shear profile.

d_volume_d_phi: jt.Float[np.ndarray, 'n_surfaces']

Derivative of plasma volume with respect to the enclosed toroidal flux phi.

well: jt.Float[np.ndarray, 'n_surfaces']

Magnetic well profile.

toroidal_current: jt.Float[np.ndarray, 'n_surfaces']

Enclosed toroidal current profile.

d_toroidal_current_d_volume: jt.Float[np.ndarray, 'n_surfaces']

Derivative of enclosed toroidal current with respect to plasma volume.

pressure: jt.Float[np.ndarray, 'n_surfaces']

Pressure profile p.

d_pressure_d_volume: jt.Float[np.ndarray, 'n_surfaces']

Derivative of pressure profile with respect to plasma volume.

DMerc: jt.Float[np.ndarray, 'n_surfaces']

Full Mercier stability criterion.

Dshear: jt.Float[np.ndarray, 'n_surfaces']

Mercier criterion contribution due to magnetic shear.

Dwell: jt.Float[np.ndarray, 'n_surfaces']

Mercier criterion contribution due to magnetic well.

Dcurr: jt.Float[np.ndarray, 'n_surfaces']

Mercier criterion contribution due to plasma currents.

Dgeod: jt.Float[np.ndarray, 'n_surfaces']

Mercier criterion contribution due to geodesic curvature.

class vmecpp.Threed1Volumetrics(**data)

Bases: BaseModelWithNumpy

Parameters:

data (Any)

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'ser_json_inf_nan': 'strings'}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

int_p: float

Total plasma pressure integrated over the plasma volume.

avg_p: float

Volume-averaged plasma pressure.

int_bpol: float

Total poloidal magnetic field energy B_phi^2/(2 mu0) integrated over the plasma volume.

avg_bpol: float

Volume-averaged poloidal magnetic field energy.

int_btor: float

Total toroidal magnetic field energy integrated over the plasma volume.

avg_btor: float

Volume-averaged toroidal magnetic field energy.

int_modb: float

Total |B| integrated over the plasma volume.

avg_modb: float

Volume-averaged |B|.

int_ekin: float

Total kinetic energy integrated over the plasma volume.

avg_ekin: float

Volume-averaged kinetic energy.

class vmecpp.MakegridParameters(**data)

Bases: BaseModelWithNumpy

Pydantic model mirroring the C++ makegrid::MakegridParameters struct.

Parameters defining the grid for external magnetic field calculations (mgrid file).

Parameters:

data (Any)

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'ser_json_inf_nan': 'strings'}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

normalize_by_currents: bool

If true, normalize the magnetic field by coil currents and windings.

assume_stellarator_symmetry: bool

If true, compute on half-period and mirror.

number_of_field_periods: int

Number of toroidal field periods.

r_grid_minimum: float

Radial coordinate of the first grid point.

r_grid_maximum: float

Radial coordinate of the last grid point.

number_of_r_grid_points: int

Number of radial grid points.

z_grid_minimum: float

Vertical coordinate of the first grid point.

z_grid_maximum: float

Vertical coordinate of the last grid point.

number_of_z_grid_points: int

Number of vertical grid points.

number_of_phi_grid_points: int

Number of toroidal grid points per field period.

static from_file(input_file)
Parameters:

input_file (str | Path)

Return type:

MakegridParameters

class vmecpp.MagneticFieldResponseTable(**data)

Bases: BaseModelWithNumpy

Pydantic model mirroring the C++ makegrid::MagneticFieldResponseTable struct.

Holds the precomputed magnetic field response on a grid, separated by coil circuit. Each field component (b_r, b_p, b_z) is a list where each element corresponds to a circuit and contains the flattened 1D array of field values on the grid for that circuit carrying unit current.

Parameters:

data (Any)

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'ser_json_inf_nan': 'strings'}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

parameters: MakegridParameters
b_r: jt.Float[np.ndarray, 'num_coils num_mgrid_cells']

Cylindrical R components of magnetic field per circuit.

b_p: jt.Float[np.ndarray, 'num_coils num_mgrid_cells']

Cylindrical Phi components of magnetic field per circuit.

b_z: jt.Float[np.ndarray, 'num_coils num_mgrid_cells']

Cylindrical Z components of magnetic field per circuit.

static from_coils_file(coils_path, makegrid_parameters)
Parameters:
Return type:

MagneticFieldResponseTable

class vmecpp.FreeBoundaryMethod(value)

Bases: str, Enum

Method for handling free-boundary conditions.

NESTOR = 'nestor'

NEumann Solver for TORoidal systems.

ONLY_COILS = 'only_coils'

Use just the coils field for the free-boundary force contribution.

This can be particularly useful for verification calculations.

Warning: This is only valid for vacuum calculations and will ignore the plasma current contribution!

BIEST = 'biest'

Boundary Integral Equation Solver for Toroidal systems.

class vmecpp.IterationStyle(value)

Bases: str, Enum

Time-step / restart control scheme for the equilibrium iteration.

VMEC_8_52 = 'vmec_8_52'

The Fortran VMEC 8.52 control (the default).

PARVMEC = 'parvmec'

The PARVMEC / VMEC2000 9.0 control.

vmecpp.set_profile(vmec_input, field, f)

Populate a line segment profile using callable f.

This allows users to set a precise pressure/iota/curent profile on flux surfaces without any precision loss by fitting the profile to e.g. experimental data first.

The callable is evaluated on all unique s values required for the multi-grid steps (full and half grids). The resulting knots and values are stored in the auxiliary arrays for the chosen profile. Therefore you should populate the multigrid ns_array resolutions before calling this function.

Parameters:

  • vmec_input (VmecInput) – The vmec input to be modified.

  • field (Literal['pressure', 'iota', 'current']) – The profile quantity to set.

  • f (Callable[[ndarray], ndarray]) – A callable taking an array of flux coordinates s and returning the value of the selected quantity (pressure/iota/current)

Return type:

VmecInput

Returns:

A modified copy of the given VmecInput.

vmecpp.iterate(vmec_input, ns=None, *, iteration_style=None, verbose=False, callback=None)

Create a single-resolution VmecModel for vmec_input and run the Python force-balance iteration on it.

iteration_style is one of "vmec_8_52", "parvmec", "robust" and overrides the input’s own setting; when None the input’s iteration_style is used (default "vmec_8_52"). Returns (model, IterationResult); model holds the converged geometry (inspect it with get_state() / the residual properties). ns defaults to the last entry of the input’s ns_array.

callback is forwarded to solve_equilibrium(): when given, it is called with an IterationState snapshot once per force iteration, for tracing the convergence and flow-control dynamics.

Parameters:
vmecpp.solve_equilibrium(model, *, style=None, verbose=False, callback=None)

Drive an initialized VmecModel to force balance.

Reproduces the C++ Vmec::SolveEquilibrium inner solve: an outer retry loop wraps the force iteration so a bad initial Jacobian triggers a magnetic axis reguess; the inner loop evaluates the forward model, tests convergence on the invariant residuals, advances the geometry with the Garabedian time step (damping from the preconditioned-residual history), and applies the time-step / restart control of the selected style. style is one of "vmec_8_52", "parvmec", "robust"; when None it defaults to model.iteration_style. model must already be initialized at a single radial resolution via VmecModel.create.

If callback is given, it is invoked once per force iteration with an IterationState snapshot of every convergence and flow-control quantity (after that iteration’s time-step / restart decision), for tracing the iteration dynamics. It is not called on iterations that immediately break out for convergence, a bad-Jacobian axis reguess, or an ijacob escalation reset.

Parameters:
vmecpp.solve_multigrid(vmec_input, *, iteration_style=None, interpolation=None, verbose=False, callback=None)

Drive the full coarse-to-fine multi-grid ramp of vmec_input from Python, mirroring the multi-grid loop of the C++ Vmec::run.

Walks ns_array in order: entries coarser than the finest resolution solved so far are skipped (the C++ ns_min rule), each remaining stage is solved to its ftol_array / niter_array tolerance with solve_equilibrium(), and the converged geometry is interpolated onto the next stage’s radial grid with VmecModel.refine_to (the C++ InterpolateToNextMultigridStep). interpolation selects the radial interpolant for those transfers: "linear" (2-point, VMEC 8.52 behavior), "cubic" (4-point Lagrange in s) or "cubic_rho" (4-point Lagrange in rho = sqrt(s)); the higher-order interpolants reduce the interpolation-added error of the stage seed down to the coarse grid’s own discretization error. None (default) keeps the linear scheme. A stage that exhausts its iteration budget without converging proceeds to the next stage, exactly as the C++ run does; a stage that hard-fails (unrecoverable bad Jacobian, the ijacob >= 75 give-up) stops the ramp.

ftol_array / niter_array entries are matched to ns_array entries by ns value (VmecModel.create / refine_to semantics); repeated ns values in ns_array are not supported because the model only refines to strictly finer grids.

Returns (model, results): the model holds the final stage’s geometry and results is the per-stage list of IterationResult.

Parameters:
class vmecpp.IterationResult(converged, failed, num_iterations, restarts, axis_reguesses, fsqr, fsqz, fsql, iteration_style='vmec_8_52', force_residual_r=<factory>, force_residual_z=<factory>, force_residual_lambda=<factory>, restart_reasons=<factory>)

Bases: object

Outcome of solve_equilibrium().

Parameters:
converged: bool
failed: bool
num_iterations: int
restarts: int
axis_reguesses: int
fsqr: float
fsqz: float
fsql: float
iteration_style: str = 'vmec_8_52'
force_residual_r: list[float]
force_residual_z: list[float]
force_residual_lambda: list[float]
restart_reasons: list[int]
class vmecpp.IterationState(iteration, iter1, iter2, fsqr, fsqz, fsql, fsq_invariant, fsqr1, fsqz1, fsql1, fsq_preconditioned, res0, res1, delt, otav, dtau, ijacob, restart_reason, n_restarts, n_reguesses, restarted, saved_backup, mhd_energy)

Bases: object

Snapshot of all convergence and flow-control quantities at one force iteration, passed to the callback of solve_equilibrium().

Captured after the time-step / restart control has run for the iteration, so delt / restarted / saved_backup reflect the decision taken this step. Use it to trace why the convergence dynamics change qualitatively (e.g. the time step collapsing, the residual leash repeatedly reverting, the damping average saturating) over the course of a run.

Parameters:
iteration: int

Global force-iteration index (1-based, monotonically increasing).

iter1: int

First iteration of the current time-step segment (reset on every restart).

iter2: int

Effective iteration index used by the forward model (iteration - bad_resets).

fsqr: float
fsqz: float
fsql: float
fsq_invariant: float

Fsqr + fsqz + fsql, the sum tested against ftolv for convergence.

fsqr1: float
fsqz1: float
fsql1: float
fsq_preconditioned: float

Fsqr1 + fsqz1 + fsql1, the sum the time-step control watches for blow-up.

res0: float

Running minimum of the preconditioned residual sum (revert leash baseline).

res1: float

Running minimum of the invariant residual sum (parvmec/robust leash).

delt: float

Current time step delt0r after this iteration’s control decision.

otav: float

Mean inverse damping time 1/tau over the last _NDAMP steps.

dtau: float

Delt0r * otav / 2, the damping applied in the Garabedian time step.

ijacob: int

Bad-Jacobian / blow-up escalation counter (give-up at >= 75).

restart_reason: int

RestartReason returned by the forward model this iteration (0 if none).

n_restarts: int

Cumulative time-step reverts so far.

n_reguesses: int

Cumulative magnetic-axis reguesses so far.

restarted: bool

True if the geometry was reverted (time step cut) this iteration.

saved_backup: bool

True if the geometry was checkpointed (both residuals improved) this step.

mhd_energy: float

Total MHD energy W of the current geometry.

class vmecpp.SolverState(iteration, multigrid_step, ns, fsqr, fsqz, fsql, ftol, delt, restart_reason, jacobian_resets, vacuum_pressure_active, mhd_energy, geometry, half_grid, curr_h)

Bases: object

The state of the solver after one force iteration, handed to the iteration_callback of run().

Parameters:
iteration: int

Iteration counter of the current multigrid stage, as printed.

multigrid_step: int

Index into ns_array of the current stage; -1 for the inserted ns = 3 stage.

ns: int

Number of flux surfaces of the current stage.

fsqr: float

Invariant force residual of R.

fsqz: float

Invariant force residual of Z.

fsql: float

Invariant force residual of lambda.

ftol: float

Tolerance the three residuals are tested against.

delt: float

Current time step.

restart_reason: int

1 no restart, 2 bad Jacobian, 3 bad progress, 4 huge initial forces; any value but 1 means the state was reverted to the last backup.

jacobian_resets: int

Jacobian resets so far in this stage.

vacuum_pressure_active: bool

Whether the vacuum pressure is part of the force balance yet.

mhd_energy: float

MHD energy of the state.

geometry: Geometry

R, Z and lambda coefficients of the state.

half_grid: HalfGridFields

The fields of the force evaluation of this iteration, which the time step that followed it has moved geometry away from, except on the iteration that converges.

curr_h: ndarray

[ns - 1] with ncurr = 1, the enclosed toroidal current that the force evaluations prescribe, in the units of half_grid.buco: each evaluation solves for chi’ so that buco equals it. Empty with ncurr = 0. The callback may change its values in place, and the force evaluations prescribe the changed current from the next iteration to the end of the multigrid step.

class vmecpp.HalfGridFields(gsqrt, bsupu, bsupv, bsubu, bsubv, weight, buco, bvco, iota, phip, vp, ntheta_even, ntheta_eff, nzeta, nfp, signgs)

Bases: object

The fields of a force evaluation on the half grid, at the angles theta_l = 2 pi l / ntheta_even, l < ntheta_eff, and zeta_k = 2 pi k / (nfp nzeta), k < nzeta.

Without lasym the poloidal points cover [0, pi], and a field on the rest of a surface follows from f(theta, zeta) = f(-theta, -zeta).

Parameters:

  • gsqrt (ndarray)

  • bsupu (ndarray)

  • bsupv (ndarray)

  • bsubu (ndarray)

  • bsubv (ndarray)

  • weight (ndarray)

  • buco (ndarray)

  • bvco (ndarray)

  • iota (ndarray)

  • phip (ndarray)

  • vp (ndarray)

  • ntheta_even (int)

  • ntheta_eff (int)

  • nzeta (int)

  • nfp (int)

  • signgs (int)

gsqrt: ndarray

the Jacobian sqrt(g), whose sign is signgs.

Type:

[ns - 1, nzeta * ntheta_eff], each row zeta-major

bsupu: ndarray

[ns - 1, nzeta * ntheta_eff] contravariant poloidal component of B.

bsupv: ndarray

[ns - 1, nzeta * ntheta_eff] contravariant toroidal component of B.

bsubu: ndarray

[ns - 1, nzeta * ntheta_eff] covariant poloidal component of B.

bsubv: ndarray

[ns - 1, nzeta * ntheta_eff] covariant toroidal component of B.

weight: ndarray

[ntheta_eff] weight of each point in an angle average, <f> = sum_{k,l} weight_l f_kl.

buco: ndarray

[ns - 1] <B_theta>, the buco of the wout file.

bvco: ndarray

[ns - 1] <B_zeta>, the bvco of the wout file.

iota: ndarray

[ns - 1] rotational transform, the iotas of the wout file.

phip: ndarray

[ns - 1] toroidal flux derivative, the phips of the wout file.

vp: ndarray

[ns - 1] signgs <sqrt(g)>, the vp of the wout file.

ntheta_even: int
ntheta_eff: int
nzeta: int
nfp: int
signgs: int
vmecpp.has_exact_force_jacobian()

Returns true if this build has the exact force Jacobian (Enzyme plugin).

This is a static build feature, cheap to query and independent of any particular equilibrium: it does not require creating a VmecModel. Use it to check ahead of time whether reverse-mode differentiation through vmecpp.autodiff is available in the current installation. A model’s has_exact_force_jacobian property additionally depends on that model’s force-balance formulation and can be false even when this function is true.

Return type:

bool

Submodules