Physics model and equations

dkx solves a radially local, steady-state drift-kinetic problem on a single flux surface. This page summarizes the governing model, the main approximations, and how the physics terms map onto the discretized operators documented in Drift-kinetic equation and system of equations, Geometry models and loading, and Numerics and algorithms. For full derivations and code-location detail, see Physics reference: the radially local drift-kinetic model and Source-code map.

Model overview

dkx evolves the non-adiabatic perturbation \(f_{s1}\) about a Maxwellian background \(f_{s0}\) on a single flux surface:

\[f_s = f_{s0} + f_{s1}.\]

In normalized variables, the linearized drift-kinetic equation (DKE) can be written as

\[\mathcal{L}_s[f_{s1}] = S_s,\]

where the operator \(\mathcal{L}_s\) includes streaming, mirror force, \(E\times B\) drifts, magnetic drifts, energy/pitch-angle drift terms, and the linearized collision operator. The source \(S_s\) contains thermodynamic drives, the inductive electric field, and the RHSMode-specific forcing used in transport matrix calculations. [1]

Background distribution and normalized variables

The code uses the normalized variables

\[x_s = \frac{v}{\sqrt{2T_s/m_s}}, \qquad \xi = \frac{v_\parallel}{v}, \qquad \mu = \frac{v_\perp^2}{2B},\]

with a background state

\[f_{s0} = n_s \left(\frac{m_s}{2\pi T_s}\right)^{3/2} \exp(-x_s^2)\]

or, when the flux-surface-varying potential is included,

\[f_{s0} = f_{sM}\exp\!\left(-\frac{Z_s e \Phi_1}{T_s}\right).\]

The normalization conventions used for hats and dimensionless drives are summarized in Normalizations and units. These conventions matter directly for the coefficients assembled by the consolidated dkx.drift_kinetic.KineticOperator and the diagnostics written by dkx.writer.

Geometry and guiding-center drifts

The guiding-center drifts can be expressed (in physical variables) as

\[\mathbf{v}_m = \frac{v_\parallel^2}{\Omega_s}\,\mathbf{b}\times(\mathbf{b}\cdot\nabla\mathbf{b}) + \frac{\mu}{\Omega_s}\,\mathbf{b}\times\nabla B,\]

and the \(E\times B\) drift as

\[\mathbf{v}_E = \frac{c}{B^2}\,\mathbf{E}\times\mathbf{B}.\]

SFINCS evaluates the geometric coefficients using Boozer-like straight-field-line coordinates (especially for geometryScheme=11/12), and the discrete operator uses those coefficients to build the drift terms in the DKE. [2]

Collision operators

SFINCS Fortran v3 supports two collision models, and dkx adds a third:

  • Pitch-angle scattering (PAS) (collisionOperator = 1): a diagonal-in-\(L\) operator used for reduced models and benchmark suites.

  • Full linearized Fokker–Planck (Landau) (collisionOperator = 0): implemented via Rosenbluth potentials and dense coupling in the speed coordinate. This is the default for high-fidelity multispecies studies.

  • Improved Sugama model operator (collisionOperator = 3): the momentum- and energy-conserving improved linearized model operator of Sugama et al. (2019), a dkx research extension beyond Fortran v3 that remains accurate into the highly collisional regime.

The linearized FP operator is the most accurate model for neoclassical transport in SFINCS and is the basis for the collision-driven preconditioners used in dkx. [1]

Constraint closure and source/sink terms

The discrete linear system is closed by auxiliary constraints. In the common linear formulation these are equivalent to enforcing

\[\left\langle \int d^3v \, f_{s1} \right\rangle = 0, \qquad \left\langle \int d^3v \, v^2 f_{s1} \right\rangle = 0,\]

with optional quasineutrality and gauge conditions when \(\Phi_1\) is solved:

\[\lambda + \sum_s Z_s \int d^3v \, (f_{s0}+f_{s1}) = 0, \qquad \langle \Phi_1 \rangle = 0.\]

These conditions remove nullspaces associated with conservation laws and determine the algebraic branch selected by the solve.

Phi1 and quasineutrality

When includePhi1 = .true., SFINCS solves for the flux-surface variation of the electrostatic potential \(\Phi_1(\theta,\zeta)\) via a quasineutrality constraint. The resulting potential modifies the kinetic equation and the collision operator through poloidal density variations. This physics is especially important for impurity transport and flows in stellarators. [3]

Transport coefficients

Transport-matrix modes (RHSMode=2/3) solve the same linearized DKE with multiple right-hand sides and postprocess the solutions into particle/heat fluxes and FSAB flows. These coefficients are the basis for neoclassical transport predictions and bootstrap current calculations in SFINCS. [1]

The main moments of interest are built from velocity-space and flux-surface integrals, schematically

\[\Gamma_s = \left\langle \int d^3v \, f_{s1}\,\mathbf{v}_{d,s}\cdot\nabla r \right\rangle,\]
\[Q_s = \left\langle \int d^3v \, \frac{m_s v^2}{2} f_{s1}\,\mathbf{v}_{d,s}\cdot\nabla r \right\rangle,\]
\[V_{\parallel,s} = \left\langle \int d^3v \, v_\parallel f_{s1} \right\rangle, \qquad j_\parallel = \sum_s Z_s e V_{\parallel,s}.\]

dkx evaluates these moments in dkx.moments (with geometry-derived scalars in dkx.magnetic_geometry).

Trajectory-model knobs

The physical model is intentionally configurable. The most important switches are:

  • magneticDriftScheme for magnetic-drift terms,

  • useDKESExBDrift for DKES-like vs full \(E\times B\) advection,

  • includeXDotTerm and the corresponding \(\dot \xi\) electric-field term,

  • collisionOperator for PAS vs full FP,

  • includePhi1 / includePhi1InKineticEquation / includePhi1InCollisionOperator for flux-surface-varying electrostatic physics.

The detailed switch-to-equation map is given in Drift-kinetic equation and system of equations.

Implementation notes

References