[spossann]

Model docs (preview): https://struphy-hub.github.io/docs-preview/pr-828/sections/subsections/models-all.html#struphy.models.DriftKineticMaxwell

The drift-kinetic counterpart of VlasovMaxwellCoupling / VlasovMaxwellQuasiNeutral (#821, merged).

What

• New propagator DriftKineticMaxwellCoupling: Maxwell's equations, the parallel acceleration of the guiding centers ($\dot v_\parallel = \frac1\varepsilon \mathbf b^* \cdot\mathbf E$) and, as an option, the $E\times B$ current in one Crank–Nicolson step at frozen positions. Gauge-invariant formulation: the guiding-center current $\int f^{\rm vol}\dot{\mathbf X}$ with $\mathbf E^* = \mathbf E - \varepsilon\mu\nabla|B_0|$, not $-\nabla\phi$. The Schur system for $\mathbf e^{\rm mid}$:
  $$\Big(\mathbb M_E + \tfrac{\Delta t^2}{4}(c_1\mathbb C^\top\mathbb M_2\mathbb C + c_2c_3\mathbb A_\parallel) + \tfrac{\Delta t}{2}c_2\mathbb S\Big)\mathbf e^{\rm mid} = \dots,\qquad \mathbb M_E = \lambda^2\mathbb M_1 + \chi\mathbb M_1^{\rm gyro}$$
  • $\mathbb A_\parallel$ (rank one per marker) is the parallel response. $\mathbb S$ (antisymmetric) is the $E\times B$ current.
  • Energy is conserved exactly in any geometry.
  • The particle matrices are assembled (dk_maxwell) or matrix-free (dk_maxwell_mf with the existing rank-1/antisymmetric MarkerOperator kernels).
  • Preconditioner: StiffnessPreconditioner. Initial projection for $\lambda = 0$, as in VlasovMaxwellCoupling: implicit energy-conserving Vlasov-Maxwell step with quasi-neutral limit; model VlasovMaxwellQuasiNeutral #821.
• Polarization term $\chi,\frac{n_0}{|B_0|^2}\partial_t\mathbf E_\perp$, using the existing M1gyro. It is needed because the drift-kinetic ordering removes the polarization drift. Without it the quasi-neutral limit $\lambda = 0$ is singular for flute modes ($k_\parallel = 0$), and with only the displacement current the shear Alfvén wave travels at the speed of light. With it, $\lambda = 0$ gives the inertial Alfvén wave. lam**2 + chi > 0 is asserted.
• PushGuidingCenterBxEstar / PushGuidingCenterParallel accept either phi (the existing path, unchanged) or a new e_field (Hcurl), with $\mathbf E^* = \mathbf E - \varepsilon\mu\nabla|B|$. The existing kernels already read the E field.
  • With e_field only the explicit algorithm is consistent: the discrete gradients conserve a Hamiltonian with a potential, and without one they would remove the work of E along the step.
  • A discrete-gradient algo therefore falls back to explicit with a warning, so the default options keep working.
• New model DriftKineticMaxwell, splitting coupling → push_para → push_bxe(e_field), and the example inertial_alfven_wave.
  • Normalized with the ion thermal velocity, as DriftKineticElectrostaticAdiabatic, using the equation parameters $\varepsilon$, $\lambda$, $\kappa$ of VlasovMaxwellCoupling: implicit energy-conserving Vlasov-Maxwell step with quasi-neutral limit; model VlasovMaxwellQuasiNeutral #821. Ampère's law becomes
    $$\varepsilon^2\lambda^2,\partial_t\mathbf E + \varepsilon^2\kappa^2,\tfrac{n_0}{|B_0|^2},\partial_t\mathbf E_\perp = \nabla\times\tilde{\mathbf B} - \varepsilon\kappa^2,\mathbf J_{gc},$$
    and the parallel force has the coefficient $1/\varepsilon$.
  • Propagator option lam (coefficient $\lambda^2$ of the displacement current; 0 is the quasi-neutral limit). Its defaults are $\varepsilon\lambda_s$, $\chi = (\varepsilon\kappa_s)^2$, $c_2 = \kappa_s^2\varepsilon$, $c_3 = 1/\varepsilon$.
  • For deuterium, 1 m, 5 T, $10^{20}$ m$^{-3}$, 10 keV: $\varepsilon^2\lambda^2 = 5\cdot10^{-6}$, $\varepsilon^2\kappa^2 = 8\cdot10^{-3}$, $\varepsilon\kappa^2 = 2.8$.
  • The constructor takes epsilon, lam and kappa as overrides. The default is BaseUnits(kBT=1.0), also in the generated params file.
• The parallel force $\frac1\varepsilon\mathbf b^*\cdot\mathbf E$ is applied only by the coupling step.
  • The model gives push_para no phi and no e_field, and asserts this.
  • doc_discretization lists which propagator applies each term.
  • A test checks that push_para leaves $v_\parallel$ unchanged under $\mathbf E \neq 0$, and changes it when given $\mathbf E$.
• The coupling propagator has no CUDA kernels yet (asserted).

Analysis note with the derivations (dispersion relation, well-posedness, energy across the splitting): notes/drift_kinetic_maxwell.md (local, not in the repo).

Energy across the splitting

The coupling step conserves energy exactly. The full splitting does so only when $\nabla|B_0| = 0$ and either $\mu = 0$ or $\tilde B_\parallel = 0$. Two pieces are missing:

• the grad-$B$ drift current, which pairs with the $E\times B$ motion of the positions (nonlinear);
• the magnetization current, which pairs with $\mu\nabla\tilde B_\parallel$ in the pushers.

Both are follow-ups.

Tests

Separate review agents went over the branch before each push. The one real finding (discrete gradients with e_field) is fixed in the second commit; the review of the normalization commit found no bugs.

• test_drift_kinetic_maxwell_coupling.py (1 and 4 ranks):

  • energy conserved to $10^{-10}$ of the exchanged energy, for $\lambda = 1, 0$, with and without $E\times B$, assembled and matrix-free, and in a Colella mapping;
  • assembled results equal matrix-free results;
  • $\mathbb A_\parallel$ is symmetric PSD and $\mathbb S$ antisymmetric;
  • the discrete equations hold;
  • the initial projection removes the period-2 oscillation;
  • dt change; invalid options.

• test_push_guiding_center_e_field.py: e_field = -grad phi reproduces the phi path exactly, for both pushers, also with the default (discrete-gradient) algo through the fallback.

• test_verif_DriftKineticMaxwell.py: frequency of the cold-slab shear mode against the roots of the dispersion determinant:

  case	measured	analytic
  $\lambda = 0$	0.3579	0.3616
  $\lambda = 0$ with $E\times B$	0.2637	0.2643
  $\lambda = 1$	0.3405	0.3428

  The deviations are marker noise (−0.35 % with 4× more markers). Energy over the whole splitting is conserved to $10^{-9}$.

• Generic model test, example test run, kernel catalog counts bumped. Regression: the generic tests of DriftKineticElectrostaticAdiabatic, GuidingCenter, ToyDrift, LinearMHDDriftkineticCC and test_current_coupling_5d.py pass.

🤖 Generated with Claude Code