Physics models

Modified

September 30, 2026

Purpose

The Tier A analytic models. At M0 they exercise the loop, the guard and the tests without any external code; they also serve as the fallbacks in toy runs. None of them is a physics claim about a real plasma. The one model meant for real scenarios, the Gaussian RF source, is at the end.

Circular equilibrium

Large-aspect-ratio circular flux surfaces with a prescribed parabolic safety factor, independent of the profiles. With \(r = a\rho\):

\[V = 2\pi^2 R_0 r^2,\quad V' = 4\pi^2 R_0 a^2 \rho,\quad G = \frac{1}{a^2},\quad q = q_0 + (q_a - q_0)\rho^2,\]

\[\frac{\partial\psi}{\partial\rho} = \frac{2\pi a^2 B_0\,\rho}{q}\ \ (\psi(0) = 0,\ \text{trapezoidal integration}), \qquad I_p = \frac{2\pi a^2 B_0}{\mu_0 R_0 q_a}.\]

For this geometry, and only for it, the toroidal-flux radius at the boundary equals \(a\) (see K-007).

Critical-gradient transport

A smooth placeholder for stiff turbulent transport (OQ-1):

\[\chi = \chi_\min + \chi_\text{stiff}\, w\, \operatorname{softplus}\!\left(\frac{R/L_T - (R/L_T)_\text{crit}}{w}\right), \qquad D_e = \frac{D}{\chi}\,\chi_e, \qquad v_e = v_\text{pinch}\,\rho .\]

Applied separately to \(\chi_e\) (from \(R/L_{T_e}\)) and \(\chi_i\) (from \(R/L_{T_i}\)). Softplus instead of a hard threshold, deliberately: the model has to be differentiable. The width \(w\) sets how stiff it is; the original M0 settings (\(\chi_\text{stiff} = 1\), \(w = 1\)) made relaxed Picard iteration oscillate, which is what motivated Newton (OQ-4). Its domain is unbounded, so it never asks for a fallback. Tested by tests/physics/test_monotonic.py (near the floor well below threshold, stiff well above it; more heating power gives higher temperatures).

Gaussian sources

Heating and particle sources with Gaussian radial shapes, normalised so the volume integral equals the configured total exactly on the discrete grid:

\[S(\rho_i) = S_\text{tot}\,\frac{\exp\!\left(-\tfrac12 ((\rho_i - \rho_\text{dep})/w)^2\right)}{\sum_j \exp\!\left(-\tfrac12((\rho_j - \rho_\text{dep})/w)^2\right) V'_j\,\Delta\rho_j}.\]

Sawtooth placeholder

The identity. A smoothed-trigger model is part of the brief’s design; runs through TORAX use TORAX’s sawtooth model where the scenario enables one (neither ITER hybrid scenario does).

Gaussian RF heating

M1 asks for an analytic Gaussian radio-frequency deposition model. It will be the trusted fallback for an RF-deposition surrogate at M3, so it is a named model rather than a one-off config. TORAX already has one (its electron-cyclotron source: a Gaussian in \(\rho\) of chosen centre and width, normalised to the requested power, all to electrons, with optional co-located current drive), so the Tokamak Toolkit’s model is a typed wrapper over it, not a second implementation:

# tkit/physics/rf.py @ ae294c9 (trimmed)
@dataclass(frozen=True)
class GaussianRF:
    p_total: float                          #: absorbed power [W]
    rho_dep: float                          #: deposition centre in rho_tor_norm [-]
    width: float                            #: Gaussian width in rho_tor_norm [-]
    current_drive_efficiency: float = 0.0   #: TORAX ECCD efficiency, 0 = heating only [-]

    def torax_source(self) -> dict[str, Any]:
        return {"P_total": self.p_total, "gaussian_location": self.rho_dep,
                "gaussian_width": self.width,
                "current_drive_efficiency": self.current_drive_efficiency}

Invalid inputs (negative power, centre outside \([0,1]\), non-positive width) are rejected at construction. Tests run the ITER hybrid scenario with 10 MW deposited at \(\rho = 0.4\):

check	result
absorbed power equals the request	10⁻⁶C-017
heating and driven current peak at \(\rho_\text{dep}\)	within one grid cell
electron temperature inside \(\rho = 0.45\) rises everywhere, and stored energy rises	yes

The last row is the brief’s “monotonic response to heating” test, on a real scenario.

Changelog

2026-09-30: first published version.