Equilibrium and kinetic profiles

Kinetic diagnostics on VEST — Thomson scattering (TS) and charge exchange spectroscopy (CES) — measure $T_e$, $n_e$, $T_i$ and $V_{tor}$ at discrete $(R, Z)$ points. To become a core_profiles IDS these point measurements must first be mapped onto a flux coordinate using an equilibrium, then fitted to a smooth 1D profile in $\rho$.

vaft.process implements that two-stage pipeline, and vaft.formula.fit_profile provides the underlying 1D fitting engine.

Pipeline

flowchart LR
    A["TS / CES .mat"] -->|machine_mapping| B["ods: thomson_scattering<br/>charge_exchange"]
    E["GEQDSK / equilibrium"] --> C
    B --> C["equilibrium_mapping_*<br/>(R,Z) → ρ"]
    C --> D["profile_fitting_*<br/>ρ → Te, ne, Ti, Vtor"]
    D --> F["core_profiles<br/>profiles_1d"]

The two stages are deliberately separate: mapping depends on the equilibrium, fitting does not. You can re-fit with a different model without recomputing the mapping.

Loading kinetic diagnostics

Raw TS/CES .mat files are turned into IDS nodes by the machine mapping layer:

import vaft
from omas import ODS
from vaft.machine_mapping.thomson_scattering import thomson_scattering
from vaft.machine_mapping.charge_exchange import charge_exchange

ods = ODS()

# thomson_scattering(ods, shotnumber, data_root=None, mat_file=None)
thomson_scattering(ods, 46051, vaft.data.data_path("legacy/46051_NeTe.mat"))

# charge_exchange(ods, shotnumber, options="ces", data_root=None, mat_file=None)
charge_exchange(ods, 47514, data_root=vaft.data.data_path("legacy/CES_47514.mat"))

data_root accepts either a directory or a path to a specific *.mat file. This populates thomson_scattering.channel[:].position.{r,z}, .t_e.data, .n_e.data (plus _error_upper siblings), and for CES charge_exchange.channel[:].ion[:].{t_i,velocity_tor}.data.

Sample files ship inside the package — see Data structures for the ODS/IDS conventions, and use vaft.data.data_path() to resolve them:

Stage 1 — equilibrium mapping

Each diagnostic channel sits at a fixed $(R, Z)$. The mapping functions interpolate $\psi(R,Z)$ from the equilibrium and return every radial coordinate the equilibrium supports, one value per channel, as a MappedPositions record:

geq = vaft.data.read_geqdsk(vaft.data.data_path("efit/g040330.00320"))

mapped = vaft.process.equilibrium_mapping_thomson_scattering(ods, geq)
mapped_ces = vaft.process.equilibrium_mapping_charge_exchange(ods, geq)

mapped.psi_norm, mapped.rho_pol_norm, mapped.rho_tor_norm   # one array each, NaN outside the LCFS
mapped.available()                                          # ('rho_tor_norm', 'rho_pol_norm', 'psi_norm')

Both accept a vaft.data.GEQDSK, an OMAS equilibrium ODS, or a legacy flux-surface mapping as geq. Channels outside the last closed flux surface are NaN in every coordinate and are dropped by the fitters. $\rho_{tor,N}$ needs the equilibrium’s $q$ profile; a legacy flux-surface mapping cannot supply it, and the record says so (mapped.rho_tor_norm_unavailable) rather than substituting $\sqrt{\psi_N}$, which is a different coordinate.

The coordinate is a choice, and the default is $\rho_{tor,N}$. Before issue #420 the mappers returned a bare $\psi_N$ array under the name “rho” and every fit was made in it. Fits are now made in the coordinate you select (coordinate=), rho_tor_norm by default; pass coordinate="psi_norm" to reproduce the old numbers. A bare array is still accepted, but only together with coordinate="psi_norm", because that is the only thing it ever meant.

The same equilibrium carried through all of its representations – the four radial coordinates side by side, a PEST grid, Miller and Fourier fits, $a/L_T$ in four gradient conventions, a prescribed island, a 3-D embedding and a camera projection, each checked against the same source, time, COCOS and flux unit – is notebooks/equilibrium_representation_reference.ipynb (#1201). Which of these is a coordinate, a projection, a representation, a derived quantity or a solver convention is set out in Equilibrium representations. Building an equilibrium is covered in Equilibrium.

Stage 2 — profile fitting

n_e_fn, T_e_fn, coeffs_ne, coeffs_te, n_e_rho, T_e_rho = \
    vaft.process.profile_fitting_thomson_scattering(
        ods,
        time_ms=320.0,
        mapped_positions=mapped,
        coordinate="rho_tor_norm",  # the default; "rho_pol_norm" and "psi_norm" are the alternatives
        Te_order=3,
        Ne_order=3,
        uncertainty_option=1,       # weight the fit by per-channel error bars
        rho_points=100,
        fitting_function_te="polynomial",
        fitting_function_ne="polynomial",
    )

Note the return order: density first, then temperature. n_e_fn and T_e_fn are FittedProfile objects: callable on any $x \in [0,1]$ of the coordinate they were fitted in, and carrying .coordinate, .method and .order so that core_profiles evaluates them on the right grid and records what produced them. n_e_rho and T_e_rho are those functions already sampled on a uniform grid of rho_points. coeffs_* is None for the gp and linear methods.

The CES counterpart mirrors it, with ion_index selecting the ion species:

Vtor_fn, Ti_fn, coeffs_vtor, coeffs_ti, Vtor_rho, Ti_rho = \
    vaft.process.profile_fitting_charge_exchange(
        ods,
        time_ms=300.0,
        mapped_positions=mapped_ces,
        Ti_order=3,
        Vtor_order=3,
        fitting_function_ti="polynomial",
        fitting_function_vtor="polynomial",
        ion_index=0,
    )

Here too the velocity function comes back before the temperature function.

Fitting methods

fitting_function_* selects the model, and every variant is dispatched to vaft.formula.fit_profile:

The (1-\rho) factor in the polynomial and exponential bases drives the fitted profile toward zero at the boundary, which is usually what you want for $T_e$ and $n_e$ but is a real assumption — gp, core_poly_edge_exp and eped_tanh do not impose it.

order only affects the polynomial-family methods; it is ignored by gp, linear and eped_tanh. eped_tanh exists for pedestal_top below rather than as a general profile fitter: it has seven parameters of its own, so passing it through the enforce_physical retry loop of the IDS fitters re-runs an identical fit at every order.

Calling the fitting engine directly

For data that is not in an IDS, use the engine itself:

y_eval, y_std_eval, fit_function, coeffs = vaft.formula.fit_profile(
    x, y, y_std,
    x_eval,
    order=3,
    uncertainty_option=1,
    fitting_function="gp",
    gp_anchor=(x_anchor, y_anchor, y_std_anchor),   # GP only
    n_restarts_optimizer=5,
)

It returns the profile on x_eval, its standard deviation (zeros for methods that carry no uncertainty estimate), a callable, and the coefficients. vaft.formula.make_fit_function(mode) builds the bare polynomial / exponential basis function if you want to fit it yourself. See Physics formulas for the rest of the formula namespace.

Pedestal top

Region reductions — core versus edge, pedestal metrics — need a boundary, and a fixed one makes results from different studies incomparable. vaft.process.profile.pedestal_top finds it from the profile instead, and records how:

from vaft.process.profile import pedestal_top

top = pedestal_top(psi_norm, p_total, quantity="p_total")
top.position   # 0.9134
top.method     # "eped_fit", or "fallback"
top.reason     # why the fallback was taken; empty when it was not
top.fit        # the FittedProfile behind it, or None

It takes arrays rather than an ODS, so it works on a profile from a kinetic-profile file just as well as on a fitted diagnostic one. quantity is required and is not defaulted: EPED defines the pedestal from total pressure, and a fit to a density or a temperature puts the top somewhere else — on one MAST discharge (45453 at 750 ms) $n_e$ gives 0.949 and $p_e$ 0.911, while $T_e$ falls back: its best tanh is 0.215 wide, which is a ramp rather than a pedestal. That is the behaviour to expect, and the reason the result carries the quantity and the method rather than just a number.

position is the tanh’s centre. The pedestal’s inner knee is inner_edge, half a width further in; which of the two a study calls “the pedestal top” differs, so both are available and neither is implied.

The fallback is PEDESTAL_FALLBACK_PSI_NORM = 0.85, and it fires when there is nothing to fit: too few points in the window, a fit that did not converge, a fitted x_ped resting on a bound, or a fitted curve that does not vary across the window by PEDESTAL_RESOLUTION_FACTOR times the residual scatter, or a tanh wider than PEDESTAL_MAX_WIDTH. That last width test matters more than it looks: the model’s own box allows a width of 0.3, and over the 171 fits of one MAST campaign 25 came back between 0.20 and 0.29 — curves that ramp across half the minor radius while reporting themselves as measured pedestals. The variation test is the interesting one — the model’s f_ped - f_sep is the $\tanh$’s asymptotic amplitude, which a wide, shallow fit through noise can make look large while the curve it actually draws is almost flat. The test is on the curve, not the parameters.

Writing core_profiles

vaft.process.core_profiles evaluates the fit callables and stores the result as a profiles_1d slice:

ods = vaft.process.core_profiles(
    ods,
    time_ms=320.0,
    mapped_positions=mapped,
    n_e_function=n_e_fn,
    T_e_function=T_e_fn,
    tol_ms=0.1,
)

The fits are evaluated on the equilibrium grid in their own coordinate and stored on grid.rho_tor_norm, grid.psi and grid.rho_pol_norm. What produced the slice is written beside it: electrons.temperature_fit.parameters records the coordinate, method and order, and core_profiles.code.parameters carries one line per slice. A slice written before #420 has no such record; that absence marks a legacy product fitted in $\psi_N$.

Thomson-only slices need an ion temperature. The Ti/Te coefficient is VEST policy, not a processing default: it lives in vest.yaml (diagnostics.core_profiles.ti_te_ratio). It is currently Ti = Te (value 1.0, σ 0.5, status assumed, so the total pressure is about twice the Thomson electron pressure), chosen on #1331 until enough CES/IDS shots exist to infer it. The earlier inference, 0.17 from 9 slices of 48224/48226/48233, is kept there as superseded_inference, as a record only. The kinetic-EFIT pipeline resolves the policy per shot with vaft.machine_mapping.core_profiles.vest_core_profiles_policy(shot) before calling core_profiles.

It writes, for the new slice index i:

• core_profiles.profiles_1d[i].time — seconds (time_ms is divided by 1000)
• core_profiles.profiles_1d[i].grid.rho_tor_norm — 100 uniform points on $[0,1]$
• .electrons.temperature (eV) and .electrons.density / .density_thermal (m⁻³)
• .electrons.temperature_fit / .density_fit — the measured channel values on their mapped $\rho$, kept alongside the fit so you can overplot data against model
• .ion[0] labelled H+, populated with the electron profiles

That last point is an explicit simplification: core_profiles copies $n_e$ and $T_e$ into the ion channel, it does not use CES data. Treat ion[0] as a placeholder unless you overwrite it yourself.

If a slice already exists within tol_ms of time_ms it is replaced rather than duplicated, so re-fitting the same time in a loop is safe.

Synthetic profiles from equilibrium pressure

When there is no Thomson data, profiles can be back-derived from the equilibrium pressure by assuming $n_e$ and $T_e$ share a shape, with $P = 2 n_e T_e e$ and $g(\rho) = \sqrt{P(\rho)/P(0)}$:

ods = vaft.omas.sample_ods()      # any ODS that holds an equilibrium

# Pin the on-axis temperature (eV) and let density follow
vaft.process.core_profiles_from_eq(ods, Te0_eV=100.0, eq_time_index=0)

# Or pin the density/temperature ratio (m^-3 per eV)
vaft.process.core_profiles_from_eq_ratio(ods, C_ne_over_Te=1.0e17, eq_time_index=0)

Both read equilibrium.time_slice[i].profiles_1d.pressure and write a core_profiles slice at the matching equilibrium time. The factor of 2 absorbs the assumption $T_i = T_e$; there is no $Z_{eff}$ or impurity modelling. These are synthetic profiles — useful to seed a code that needs some kinetic input, not a measurement.

Analytic L-mode, H-mode and ITB states

The opposite direction — prescribe the kinetic profiles, derive the pressure — is a small preset layer (#1045). Each channel is a generalized-parabolic core plus optional Groebner tanh steps (the #552 kernels): an edge pedestal whose top value holds exactly at the knee, and an internal transport barrier with its own position, full width and step amplitude per channel. Axis and separatrix values hold exactly; n_i follows from quasi-neutrality with a uniform z_eff, and p_e, p_i, p and dp/dψ_N are derived, never specified. The coordinate is psi_norm; units are m⁻³, eV and Pa.

from vaft.process.equilibrium import solovev_example
from vaft.process.profile import analytic_hmode_itb_state, project_plasma_state

state = analytic_hmode_itb_state(te_axis=400.0, te_ped=80.0, te_sep=10.0,
                                 pedestal_position=0.93, pedestal_width=0.05,
                                 itb_position={"n_e": 0.3, "T_e": 0.3, "T_i": 0.45})
geometry = solovev_example("limited")
p_rz = project_plasma_state(state, geometry, "p_total")   # NaN outside the LCFS

The presets are analytic_lmode_state, analytic_hmode_state, analytic_itb_state and analytic_hmode_itb_state; compose_analytic_profile and compose_plasma_state build any other combination, and vaft.plot.plot_plasma_state_projection draws a projected quantity. A state is not an equilibrium: projecting it onto a geometry leaves that geometry solved for its own sources. See notebooks/analytic_plasma_state_presets.ipynb.

Synthetic kinetic profiles from an equilibrium

generate_synthetic_kinetic_profiles (#122) generalizes both directions above: it takes an equilibrium (GEQDSK, ODS or EquilibriumData) and a SyntheticKineticSpec of declared assumptions, and returns n_e, T_e, T_i, the main-ion and impurity densities and the pressures on the equilibrium’s psi_norm grid, with rho_tor_norm derived from its q. The fidelity ladder is explicit in the records (vaft.data): Level 0 SqrtPressureSplit (the core_profiles_from_eq closure above, bit for bit), Level 1 ProfileSpec with an AnalyticProfile (the #1045 composition) or a TabulatedProfile, Level 2 a ScalarTarget (axis, separatrix, line or volume average, Greenwald fraction), a peaking factor, T_i/T_e or p_e/p, and Level 3 a GradientProfile of prescribed a/L. The composition is a hydrogenic main ion plus at most one impurity at uniform Z_eff.

from vaft.data import (Composition, ProfileSpec, ScalarTarget, SyntheticKineticSpec,
                       TemperatureAssumption)
from vaft.data.resources import sample_geqdsk
from vaft.process.profile import (compose_analytic_profile, generate_synthetic_kinetic_profiles,
                                  write_synthetic_core_profiles)

shape = compose_analytic_profile("n_e", axis_value=1.0, separatrix_value=0.2, core_beta=1.5)
spec = SyntheticKineticSpec(
    # peaking is solved through the core exponent (beta >= 1), so it cannot go below the beta = 1
    # value of this shape (about 1.74 here); a lower request comes back constraint_not_reached.
    n_e=ProfileSpec(shape, ScalarTarget("greenwald_fraction", 0.3), peaking_factor=1.8),
    temperature=TemperatureAssumption(ti_over_te=0.5),
    composition=Composition("D", "C", z_eff=2.0),
    pressure_constraint="equilibrium", closure="temperature",   # solve T_e so p_kin = p_eq
)
result = generate_synthetic_kinetic_profiles(sample_geqdsk(), spec, time=0.319)
result.status, [(t.name, t.requested, t.achieved) for t in result.targets]
ods = write_synthetic_core_profiles(result)                      # refused unless status == "success"

pressure_constraint="equilibrium" makes p_kin = p_eq locally and names the solved variable (closure="temperature", "density" or "sqrt_split"); "thermal_energy" matches only W = 3/2 ∫ p dV and reports the local mismatch; "kinetic" keeps the profiles and reports it. Every target is recomputed from the final arrays, a contradictory request (a prescribed T_e with a closure that solves T_e) raises SyntheticProfileError with a status, and an incompatible one returns pressure_closure_failed. The profiles are assumptions, never measured or transport-predicted, and the slice says so. See notebooks/synthetic_kinetic_profiles_from_equilibrium.ipynb.

A self-consistent equilibrium and kinetic state

vaft.code.chease_kinetic_iteration.build_consistent_state (#123) makes the equilibrium and the synthetic profiles agree under a declared EquilibriumKineticSpec. With closure="equilibrium_pressure" the equilibrium pressure is authoritative: #122 decomposes it once, and CHEASE is not called. With closure="kinetic_pressure" the kinetic pressure is authoritative. CHEASE re-solves the fixed-boundary equilibrium with PRES/PPRIME built from p_kin (a PCHIP interpolant and its derivative) and FF' from an explicit CurrentPolicy: the initial FF' shape or an analytic g, with I_p or q95 held. The profiles are then regenerated on the new equilibrium from the same spec, and the loop repeats. Every state is kept with its residuals and validation, and the status names the outcome: converged, max_iterations, diverged, oscillating, stagnated, or a failure. The result is an assumption-driven self-consistent state, not a transport prediction or a reconstruction.

from vaft.code.chease_kinetic_iteration import build_consistent_state
from vaft.data import (EquilibriumKineticSpec, ProfileSpec, SyntheticKineticSpec,
                       TemperatureAssumption)
from vaft.data.resources import sample_geqdsk
from vaft.process.profile import compose_analytic_profile

n_e = ProfileSpec(compose_analytic_profile("n_e", axis_value=1e19, separatrix_value=2e18))
decompose = EquilibriumKineticSpec(
    SyntheticKineticSpec(n_e=n_e, temperature=TemperatureAssumption(ti_over_te=0.5),
                         pressure_constraint="equilibrium", closure="temperature"),
    closure="equilibrium_pressure")                              # no CHEASE solve
state = build_consistent_state(sample_geqdsk(), decompose, time=0.319)
state.status, state.initial.metrics["closure_max_relative"], len(state.iterations)

See notebooks/self_consistent_equilibrium_kinetic_iteration.ipynb for the kinetic_pressure loop on 39915.

Plotting

Each stage has a canonical plot, reached through the vaft.omas.plot_* adapters. They draw what the ODS already holds — fitting is a vaft.process step, not a plotting option:

vaft.omas.plot_thomson_scattering_geometry_poloidal(ods)            # channel positions in the poloidal plane
vaft.omas.plot_thomson_scattering_time_electron_temperature(ods)    # per-channel Te history
vaft.omas.plot_thomson_scattering_time_electron_density(ods)        # per-channel ne history
vaft.omas.plot_thomson_scattering_profile_electron_temperature(ods) # measured Te versus position
vaft.omas.plot_electron_temperature_profile(ods)                    # the fitted core_profiles Te
vaft.omas.plot_core_profiles_time_volume_averaged(ods)
vaft.omas.plot_equilibrium_profile_pressure(ods)                    # consistency check

plot_equilibrium_profile_pressure is the useful sanity check after a fit: if the pressure implied by the kinetic profiles disagrees badly with the equilibrium pressure, the fit or the mapping is wrong.

Measurement and fit are separate plots now. plot_thomson_scattering_profile_electron_temperature shows the channels; plot_electron_temperature_profile shows whatever core_profiles holds — a fit from vaft.process.profile_fitting_thomson_scattering, or the synthetic profiles written by the core_profiles_from_eq* helpers above. It takes coordinate= — rho_tor_norm by default, or psi_norm. (The equilibrium profiles accept more, including r_major; a core-profile plot does not.)

For CES, map the channels onto the equilibrium and fit in vaft.process, then plot:

mapped = vaft.process.equilibrium_mapping_charge_exchange(ods, geq)
fit = vaft.process.profile_fitting_charge_exchange(
    ods, time_ms=300.0, mapped_positions=mapped, ion_index=0,
    fitting_function_ti="polynomial", Ti_order=3,
)

vaft.omas.plot_charge_exchange_profile_ion_temperature(ods)   # measured Ti versus position
vaft.omas.plot_charge_exchange_time_ion_temperature(ods)      # per-channel Ti history

Normalized gradients

plot_electron_temperature_profile_gradient, plot_electron_density_profile_gradient and plot_ion_temperature_profile_gradient draw $-L\,d\ln f/dx_g$ of the stored core_profiles slice (#551). The builder calls vaft.process.profile_gradients.profile_gradient once, and the renderer only draws the result. Four options keep the choices apart:

The default is therefore $a/L$, and the axis label is generated from the resolved record: $a/L_{T_i}$, $R_0/L_{T_i}$ for R_major_axis, or $-\partial\ln T_i/\partial r$ [m$^{-1}$] for "none".

import vaft

ods = vaft.omas.load(vaft.data.sample(48224))
vaft.omas.plot_ion_temperature_profile_gradient(ods, coordinate="rho_tor_norm", convention="tglf")
model = vaft.omas.extract_ion_temperature_profile_gradient(ods, reference_length="R_major_axis")
model.y_label                                  # '$R_0/L_{T_i}$'
model.metadata["mathematical_definition"]      # '-R_0 * d(log(T_i)) / d(r_minor)'
model.to_xarray().attrs["metadata"]            # the same record, as JSON text

The profile is paired with the equilibrium slice at its own time, never by index. Its grid is read from grid.psi (normalized by that slice’s axis and boundary flux), else grid.rho_pol_norm, else grid.rho_tor_norm, and the record names the one used. A stored grid.rho_tor_norm is used only when it equals the slice’s own profiles_1d.rho_tor_norm, is not the $\sqrt{\psi_N}$ proxy, and lies within RHO_TOR_GRID_TOLERANCE (0.05) of the map’s $q$-integrated $\rho_{tor,N}$; the check is recorded in profile_grid_check. Samples with $f \le 0$ or a non-finite value at either end of the profile are left out and counted by side and cause in excluded_points and in the title. Everything else is refused with its reason:

• an equilibrium that is not at the profile’s time;
• a non-positive value inside the profile;
• a point outside the map’s support;
• gs2, gkw, gene and L_ref, which need run settings a plot option cannot carry. For these, call profile_gradient with metadata=;
• a flux-label gradient_coordinate with a reference length (pass reference_length="none");
• an explicit reference_length=None, which profile_gradient reads as "none".

Choosing a convention in the interactive controls stops the gradient_coordinate and reference_length controls from being sent, since the preset resolves both.

Kinetic-profile files

vaft.data.kinetic_profiles is the container the kinetic-profile file formats read into and write from — GPEC’s .kin today, the Osborne pfile and MARS PROF*.IN to follow:

from vaft.data import read_kin, write_kin, normalize_psi

profiles = read_kin("g045453.00750.kin")
profiles.psi_norm[0]              # 0.00494621873 — exactly what the file said
profiles.available()              # ('n_e', 'n_i', 'T_e', 'T_i', 'omega_exb')
profiles.normalization.method     # "as_read"
write_kin(profiles, "again.kin")  # byte-identical to this file; see below

A file VAFT wrote round-trips byte for byte, and so do the files the MAST-U workflow produces — 59 of the 65 .kin files in the reference tree, the named example among them. The other six are not VAFT’s: GPEC’s bundled DIII-D example right-aligns its columns, so a negative rotation eats a separator space, and the MARS-input file uses no leading indent at all. Both read correctly; they are simply written back in this writer’s layout.

The container fixes one unit set — densities in m⁻³, temperatures in eV, angular frequencies in rad/s, pressures in Pa — because those are what GPEC’s own reader consumes. A reader converts into them; nothing carries free-text units around.

The radial coordinate is never rescaled on read. Making it span exactly [0, 1] is normalize_psi(), an explicit operation whose result records that it happened. It has no default method=: passing axis=/edge= and getting a min–max stretch instead would be the same silent rescale, only now with a provenance record vouching for it. This is not fastidiousness: GPEC’s read_kin re-splines a .kin onto a uniform 101-point [0, 1] grid with extrapolation (nkin = 100 intervals), so it expects a truncated edge and handles it. A real MAST-U file spans ψ_N = 0.00495 → 1.0, so stretching it moves every interior point — and once the stretched values are written back, permanently.

Toroidal rotation and the E×B frequency are separate fields. omega_tor and omega_exb are never merged, and write_kin refuses a profile set carrying only the former: the .kin rotation column is ω_E, and a toroidal rotation written there is wrong in a way nothing downstream detects. write_kin also refuses to write a file GPEC would read as something other than what it says: omega_exb holding any zero (GPEC substitutes 1e-9 element by element — pass allow_zero_rotation=True for a set converted from an all-zero PROFROT.IN, which sample/output/converted_from_transp.kin is), a value that is not finite, or a radial coordinate that does not increase. Line endings are LF, so a CRLF file does not round-trip byte-identically.

Osborne pfiles

A pfile is read in two steps, because it carries more than a profile set does: per-section units, a derivative column, and an N Z A of ION SPECIES block.

from vaft.data import read_pfile, write_pfile, kinetic_profiles_from_pfile

pf = read_pfile("p045453.00750")
pf.keys()                        # the 22 sections, in file order
pf.unit("te")                    # 'KeV' — what the file declares, not what we assume
pf.section("omgeb").derivative   # the file's own third column, kept

write_pfile(pf, "again")         # byte-identical: see below

profiles = kinetic_profiles_from_pfile(pf)   # now in m^-3, eV, rad/s, Pa, V/m

PFile is the file as written and converts nothing; kinetic_profiles_from_pfile applies the one unit ladder. Writing puts the sections back into the format’s own order, so a file that was already canonical — as all 57 reference files are — comes back byte for byte, and one that was not is normalised rather than reproduced. Line endings are LF. The conversion is deliberately lossy in one direction: KineticProfiles holds profiles, so the derivative columns and the per-section units stay on the PFile, which is where a future pfile writer has to take them from. Keeping them apart means “this is byte-for-byte the file we read” and “these are the right units” can fail independently.

The derivative column is data, not something to recompute. Writing it back as read reproduces all 57 reference files byte for byte; recomputing it with np.gradient reproduces none of them — at most 8 of a file’s 22 sections agree with a recomputation to one part in a million, and no section agrees across every file. write_pfile therefore preserves it and computes one only for a section built in memory, or when you pass recompute_derivatives=True because you changed the values.

The rotation family is ten sections whose names differ by two letters and whose meanings do not. Three get fields of their own, and the mapping is by exact section name:

The other seven — omgvb, omgpp, ommvb, ommpp, omevb, omepp, omghb — keep their pfile names in extras, in the file’s own units, with those units recorded in provenance, as do kpol, vtor1 and vpol1. ptot maps to p_total as the file’s own total, fast ions included, and pb to p_fast; T_z has no pfile source.

A section on a different ψ column from the rest is refused rather than quietly reprojected — all 22 sections share one coordinate bit-for-bit in every reference file — and so is a short or ragged section. An unrecognised unit raises on a mapped section, where guessing a factor would be a silent factor of a million, but not on one kept in extras.

MARS PROF*.IN

MARS reads its kinetic profiles as a deck of two-column ASCII files, one quantity each.

from vaft.data import read_mars_profiles, write_mars_profiles

profiles = read_mars_profiles("mars_input/")      # the whole deck
profiles.omega_tor                                 # PROFROT.IN, rad/s
profiles.omega_exb                                 # PROFWE.IN, rad/s
profiles.normalization.method                      # "mars_s_squared" — see below
write_mars_profiles(profiles, "out/")

The header’s second field names the abscissa, and it is not ψ_N. MARS branches on it in every one of its profile readers and stops on anything else: key 1 means the column is s, the square root of the normalised poloidal flux, and key 2 means the toroidal equivalent (marsq.f:16269, “RAD = SQRT(TOROIDAL FLUX)”). So a key-1 deck’s column is s, and psi_norm is s² — reading it straight through puts every point at the square root of where it belongs, so 0.5 in the file becomes ψ_N = 0.25. The conversion is recorded as normalization.method == "mars_s_squared", the writer converts back (sqrt(psi_norm), header key 1), and a key-2 deck is refused: reaching poloidal flux from a toroidal-flux coordinate needs an equilibrium, which a file-format reader has not got.

The two rotation files are two different quantities, and this is where C-27 finally closes. MARS’s own source settles which is which: PROFROT.IN is the bulk toroidal fluid rotation ω_φ, read under NPROFR, and PROFWE.IN is the toroidal E×B rotation ω_E, read under NPROFWE into a different array. MARS states the relation between them in its own analytic branch — ROTWE = ROT - OMEGAI*, i.e. they differ by the ion diamagnetic frequency. So there is no precedence rule here: a deck carrying both reads as two fields, and write_mars_profiles will not fill PROFROT.IN from omega_exb, which is what produced every committed deck’s byte-identical pair.

PROFDEN.IN fills n_e and nothing else — MARS carries one density, so equating the ion density with it is a modelling choice a caller makes. Any other PROF*.IN is kept in extras under its stem.

Nothing is scaled, which is a finding rather than an assumption: “a MARS input is Alfvén-normalised” is the obvious wrong guess. The committed profiles run 1.4e4–1.0e5 rad/s against an Alfvén frequency of order 1.4e6 for MAST, and MARS converts the file itself when its NEXPV flag says the file carries absolute values — multiplying the on-axis value by the Alfvén time, which is only dimensionally sensible if the file is in rad/s. Note the corollary: under MARS’s default NEXPV = 0 only the shape of each profile is used, the scale coming from a namelist amplitude, so whether a deck’s numbers reach a run at all is a property of its RUN.IN rather than of these files.

Files whose ψ columns disagree, in value or in length, are refused rather than reprojected onto one another, and so is a header row count that disagrees with the rows beneath it.

Exporting

To hand fitted electron profiles to an external code:

vaft.process.export_electron_profile_txt(
    n_e_fn, T_e_fn, coeffs_ne, coeffs_te,
    rho_points=100,
    filename="electron_profiles.txt",
)

The file is CSV with the header psi_N, T_e [eV], n_e [m-3].

Full example

import numpy as np
import vaft
from omas import ODS
from vaft.machine_mapping.dataset_description import dataset_description
from vaft.machine_mapping.thomson_scattering import thomson_scattering

shot = 40330

geq = vaft.data.read_geqdsk(vaft.data.data_path("efit/g040330.00320"))
ods = geq.to_omas()
dataset_description(
    ods, source=shot, options={"source_type": "shot"},
)
thomson_scattering(
    ods, 46051, vaft.data.data_path("legacy/46051_NeTe.mat"),
)

mapped = vaft.process.equilibrium_mapping_thomson_scattering(ods, geq)

for t_s in np.asarray(ods["thomson_scattering.time"], dtype=float):
    time_ms = float(t_s) * 1e3
    n_e_fn, T_e_fn, *_ = vaft.process.profile_fitting_thomson_scattering(
        ods, time_ms, mapped,
        Te_order=2, Ne_order=2,
        fitting_function_te="polynomial",
        fitting_function_ne="polynomial",
    )
    ods = vaft.process.core_profiles(ods, time_ms, mapped, n_e_fn, T_e_fn)

vaft.omas.plot_thomson_scattering_profile_electron_temperature(ods)

Note that thomson_scattering.time is in seconds while the fitting and core_profiles functions take time_ms in milliseconds — the conversion above is not optional.

Choosing a fitting method

VEST Thomson systems have few channels, so the fit is under-constrained and the method choice matters more than it would on a large device:

• Start with polynomial at order=2. Higher orders oscillate between sparse channels.
• Use gp when you want an uncertainty band rather than a point estimate — it is the only method that returns a meaningful y_std_eval.
• Use core_poly_edge_exp when the edge gradient is the quantity of interest, since the polynomial bases force the profile to zero at $\rho = 1$ regardless of the data.
• linear is a diagnostic aid, not a physics result: it interpolates the channels exactly and will happily reproduce measurement noise.

Always keep uncertainty_option=1 if error bars exist. With it disabled every channel is weighted equally, including ones the diagnostic itself flagged as unreliable.

Equilibrium reconstruction with EFIT

EFIT’s process-level result.ok is retained for compatibility. For scientific use, inspect result.usable and result.slice_statuses: each time slice separately reports runtime, output, numerical, and physical status, with stable codes in vaft.code.EFIT_FAILURE_CODES. A status can be serialized through to_dict() and from_dict().

Typed EFITConfig, EFITProfileConfig, and EFITNumericsConfig make routine k-file settings explicit. prepare_efit_inputs(ods, config) records the resolved settings, stable hash, VAFT version, provenance, and k-file checksums in efit_configuration.json. For deterministic parameter scans that do not run the EFIT binary, vaft.code.efit_parameter_grid() accepts dotted paths such as profile.kppcur and constraints.group_weights.bpol_probe.

The profiles above are mapped on an equilibrium, and on VEST that equilibrium is an EFIT reconstruction: vaft.code.efit writes the k-files, runs the code and collects its output. The per-shot stages live in Automated pipelines and the functions in the code API. Two parts of that machinery belong on this page: the named configurations a run selects, and the poloidal field derived from the flux map it returns.

Named configurations (presets)

vaft.code.efit.PRESETS names complete EFIT configurations: the scientific configuration written into the k-file, sigma floor included. Since 2026-10-01 the defaults of EFITScientificConfig are the #891 working setting, so a run that names no configuration writes it (vaft.code.efit.DEFAULT_PRESET == "statistical_891").

generate_kfile applies the configured floor to a copy of the constraints, so the product itself is never changed. prepare_constraints shows what the floor does, per slice and family:

from vaft.code.efit import efit_preset

preset = efit_preset("statistical_891")
ods = vaft.omas.sample_ods()   # shot 39915, whose equilibrium carries the EFIT constraints
constraints, floor_changes = preset.prepare_constraints(ods)   # a floored copy of the constraints ODS
floor_changes[0]   # {'slice': 0, 'family': 'bpol_probe', 'floor': ..., 'raised': ..., 'fitted': ...}

Writing the k-files needs no configuration for the default; the legacy one is named:

from vaft.code.efit import generate_kfile

generate_kfile(ods, 39915, save_dir="efit-run")                                       # statistical_891
generate_kfile(ods, 39915, save_dir="efit-legacy", config=efit_preset("routine").scientific)

routine_profile_config(), routine_numerics_config(), routine_constraint_config() and routine_scientific_config() in vaft.code.efit.config build the legacy pieces one at a time. The legacy configuration is selected only by name: a positional npprime/nffprime on generate_kfile without a config (deprecated; it warns once) and EFITConfig.npprime/nffprime both mean the same thing, a basis override on the configuration otherwise in force – the default when nothing else is named. To run the legacy configuration with another basis, pass routine_scientific_config(profile=routine_profile_config(kppcur=..., kffcur=...)); to run it through prepare_efit_inputs, pass its profile, initialization, numerics and constraints.

Selecting a preset in the pipelines.

• Pipeline 1: efit.preset in config.yaml is empty for the default and routine for the legacy configuration. The k-file stage always writes efit_preset.json beside its manifest, and the EFIT product carries that record under code.parameters (efit_collection.efit_preset). vaft.database.summary reports it as efit_configuration (name@sha12, or unrecorded for products written before the record existed).
• Pipeline 2: kinetic.efit_preset builds the kinetic lineages’ base magnetic k-file the same way (KineticEFITConfig.efit_preset).

A routine run writes the same product paths as a default one, so give it its own base_dir.

The default allows up to 514 iterations (routine: 100). A slice that never converges runs to that cap, and pipeline 1 runs all of a shot’s slices in one EFIT call, so efit.timeout is 3600 s: 600 s timed out 26 of 33 shots in the #1331 Tier A campaign.

statistical_891 was calibrated on 39915’s flat-top. The weight-study README (workflow/efit_uncertainty_calibration) records how.

Poloidal field on the boundary and the flux convention

The poloidal field on the LCFS, which shafranov_integrals turns into $S_1, S_2, S_3$ for $\beta_p$ and $l_i$, is read off the flux map by bicubic spline. The map’s flux convention is an argument, not an assumption: an ODS stores $\psi$ in Wb (COCOS 11), while the function’s default is Wb/rad, a factor $2\pi$ in the field.

from vaft.process.equilibrium import poloidal_field_at_boundary

ods = vaft.omas.sample_ods()   # shot 39915: psi stored in Wb, and the record says so
slice_0 = ods["equilibrium.time_slice.0"]
R_grid_1d = slice_0["profiles_2d.0.grid.dim1"]
Z_grid_1d = slice_0["profiles_2d.0.grid.dim2"]
psi_grid = slice_0["profiles_2d.0.psi"]
R_bdry, Z_bdry = slice_0["boundary.outline.r"], slice_0["boundary.outline.z"]

B_p_bdry, B_R_bdry, B_Z_bdry = poloidal_field_at_boundary(
    R_grid_1d, Z_grid_1d, psi_grid, R_bdry, Z_bdry,
    cocos=vaft.omas.ods_cocos(ods),  # the flux convention the record declares: 11 for IMAS psi in Wb; omitting it assumes Wb/rad
)

The rest of the chain – shafranov_integrals, efit_virial_volume_integrals, calculate_diamagnetism – is documented in the process reference.

References

• Notebook: profile_fitting_using_equilibrium_and_kinetic_diagnostics.ipynb
• Source: vaft/process/profile.py
• Source: vaft/plot/profile.py
• Source: vaft/formula/utils.py
• Batch pipeline: workflow/automatic_pipeline_2_corrective_data_update/legacy/update_thomson_scattering_and_core_profile.py
• Related: Equilibrium · Data structures · Physics formulas · Examples