Drift-kinetic equation and system of equations

dkx starts from the radially local, steady drift-kinetic equation on one flux surface. The unknown is the non-adiabatic part of the species distribution, f_s1; outputs such as particle flux, heat flux, parallel flow, bootstrap-current proxies, NTV, and transport-matrix coefficients are velocity-space and flux-surface moments of this solved perturbation.

The continuous equation is treated by expanding pitch-angle dependence in a Legendre representation, discretizing speed on the x grid, discretizing theta/zeta on a periodic flux-surface grid, and applying finite-difference or spectral-style periodic stencils to the streaming, drift, and mirror terms. The collision block is either the pitch-angle-scattering model or the full linearized Fokker-Planck/Rosenbluth operator, depending on collisionOperator. Only after this discretization does the problem become the block linear/nonlinear algebraic system described below.

In its most general supported configuration, dkx solves a coupled system consisting of:

  1. a drift-kinetic equation (DKE) for each kinetic species,

  2. an optional quasineutrality equation for the flux-surface variation of the electrostatic potential \(\Phi_1(\theta,\zeta)\),

  3. auxiliary constraints that remove nullspaces and enforce moment conditions.

This page summarizes the block structure and the equations most relevant to the linear system that dkx assembles and solves. For physics background and literature references, see Physics model and equations and Physics reference: the radially local drift-kinetic model.

For full upstream context and derivations, see the SFINCS v3 technical notes and user manual in the upstream SFINCS project repository, cited from Upstream SFINCS sources and primary literature.

Unknown ordering

The operator uses a global ordering for the state vector and the rows/columns of the master system that follows the mature SFINCS-style block layout because that structure is physically natural and convenient for diagnostics, comparison, and testing.

For the common readExternalPhi1 = .false. case:

  • The F-block (distribution function) is ordered first, with indices nested as:

    species → x → xi/Legendre mode → theta → zeta.

  • If includePhi1 = .true., the Phi1 block (labeled BLOCK_QN in Fortran) contributes \(N_\theta N_\zeta\) additional unknowns corresponding to \(\Phi_1(\theta,\zeta)\) on the grid.

  • A final scalar unknown lambda enforces the constraint \(\langle \Phi_1 \rangle = 0\).

  • Constraint-scheme-dependent source unknowns (and their corresponding constraint rows) are appended last.

Linearized drift-kinetic equation (RHSMode=1)

SFINCS v3 solves a linearized, steady-state drift-kinetic equation for the non-adiabatic part of the distribution function on a single flux surface. We write

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

where \(f_{s0}\) is the Maxwellian (optionally modified by \(\Phi_1\)) and \(f_{s1}\) is the unknown. In normalized v3 form, the matrix-free operator assembled by dkx corresponds to

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

with

\[\mathcal{L}_s = \underbrace{v_\parallel \mathbf{b}\cdot\nabla}_{\text{streaming}} +\underbrace{\mathcal{M}}_{\text{mirror}} +\underbrace{\mathbf{v}_{E}\cdot\nabla}_{E\times B} +\underbrace{\mathbf{v}_{m}\cdot\nabla}_{\text{magnetic drifts}} +\underbrace{\dot{x}\,\partial_x + \dot{\xi}\,\partial_\xi}_{\text{energy / pitch-angle drifts}} -\underbrace{\sum_b \mathcal{C}^{\mathrm{lin}}_{sb}}_{\text{collisions}},\]

and a source term \(S_s\) that includes the background thermodynamic drives (\(\partial_\psi n_s\), \(\partial_\psi T_s\)), the inductive electric field (\(E_\parallel\)), and (in transport-matrix modes) the v3 whichRHS overwrites.

The explicit expressions for the ExB and magnetic-drift prefactors are summarized in Method overview, with longer-form derivations in Physics reference: the radially local drift-kinetic model.

Term-to-input switch mapping

The equation above is the full v3 linear operator. Each term is enabled (or modified) by specific namelist parameters:

  • Parallel streaming + mirror: always included (in Legendre basis, couples \(L\leftrightarrow L\pm 1\)).

  • ExB drift terms (\(\mathbf{v}_E\cdot\nabla\) in \(\theta,\zeta\)): controlled by useDKESExBDrift (DKES vs default v3 coefficient form) and by the radial electric field input (dPhiHatdpsiHat).

  • Magnetic drift terms (\(\mathbf{v}_m\cdot\nabla\) + associated \(\partial_\xi\) term): enabled when magneticDriftScheme > 0 (schemes 1–9 are fully supported in dkx); upwinding is controlled by magneticDriftDerivativeScheme.

  • Energy derivative term (\(\dot{x}\,\partial_x\)): enabled with includeXDotTerm = .true..

  • Pitch-angle drift term (\(\dot{\xi}\,\partial_\xi\)): included in v3 with the ExB and magnetic drift models; its coefficient matches populateMatrix.F90 for the supported drift schemes.

  • Collisions: choose with collisionOperator:

    • collisionOperator = 1 → pitch-angle scattering (PAS), diagonal in \(L\).

    • collisionOperator = 0 → full linearized Fokker–Planck operator (Rosenbluth form), dense in \(x\) and coupled across species.

    Additional collision options:

    • includeTemperatureEquilibrationTerm: adds the \(C_{ab}[f_{aM}, f_{bM}]\) term to the RHS.

    • includePhi1InCollisionOperator (requires includePhi1InKineticEquation): modifies collision coefficients with the poloidal density factor \(\exp(-Z_s \alpha \Phi_1/T_s)\).

  • Phi1 coupling in the kinetic equation: includePhi1 = .true. enables the QN block; includePhi1InKineticEquation = .true. inserts the Phi1-dependent terms in the DKE. If readExternalPhi1 = .true., \(\Phi_1\) is read from file and the kinetic system is solved linearly with fixed \(\Phi_1\).

  • Quasineutrality + adiabatic response: includePhi1 + quasineutralityOption select the QN equation, and withAdiabatic = .true. adds adiabatic species contributions.

Transport-matrix modes (RHSMode=2/3)

For RHSMode=2 (transport matrix) and RHSMode=3 (monoenergetic transport), v3 loops over whichRHS and overwrites selected drives internally before forming the RHS:

  • dnHatdpsiHat / dTHatdpsiHat (thermodynamic gradients),

  • EParallelHat (inductive field),

  • and, for RHSMode=3, the single-point x=1 grid used to compute monoenergetic coefficients.

On the canonical stack this behavior is KineticOperator._with_rhs_settings, which overwrites the whichRHS drive columns before KineticOperator.rhs builds each RHS, so the transport-mode right-hand side matches v3 evaluateResidual(f=0) exactly.

Quasineutrality and Phi1 constraint

The coupled system includes a quasineutrality condition and the constraint \(\langle \Phi_1 \rangle = 0\):

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

In the fully nonlinear v3 configuration, the quasineutrality equation contains additional nonlinear dependence through \(f_{s0}(\Phi_1)\) and (depending on options) adiabatic responses.

Implementation note

Not every optional Phi1 coupling that appears in the extended literature is active in every public workflow. The supported scope is documented in Inputs (namelist) reference, Outputs (HDF5, NetCDF4, and NPZ), and Validation against reference implementations. On the canonical stack the blocks above are built by the consolidated dkx.drift_kinetic.KineticOperator (streaming/mirror, ExB, the \(E_r\) terms, the quasineutrality rows, and the <Phi1>=0 Lagrange row), by dkx.collisions (PAS and Fokker–Planck blocks), and by dkx.phi1 (the nonlinear \(\Phi_1\) Newton solve). See the equation-to-code map in Physics reference: the radially local drift-kinetic model.