[spossann]

Summary

New model VlasovPoissonElectrostaticAdiabatic

The full-orbit (6D) counterpart of DriftKineticElectrostaticAdiabatic: Vlasov ions in the electrostatic potential and a static B0, electrons through the linearized Boltzmann response.

df/dt + v . grad f + (1/eps) (-grad phi + v x B0) . df/dv = 0
-(eps/alpha^2) Delta phi + n0/(Z eps T0) (phi - <phi>) = int f d^3v - n0

• New propagator AdiabaticPoissonSolve (subclass of ImplicitDiffusion):
  • quasi_neutral=True (default) drops the vacuum term, so the solve is a weighted mass-matrix solve;
  • which_geometry adds the optional flux-surface average phi - <phi>, as in the drift-kinetic model;
  • stab_mat="M0ad_withT" (weight n0/T0).
• ImplicitDiffusion now allows stab_mat="M0ad_withT" and diffusion_mat=None (no diffusion term), with option checks.
• Propagator order push_eta, poisson, push_v_E, push_vxb: with Strang splitting both E-kicks use the potential of the current marker positions (drift-kick-kick-drift). This is the same reasoning as DriftKineticElectrostaticAdiabatic: Poisson solve between the two guiding-center pushers #799 for the drift-kinetic model.
• Tests:
  • manufactured solutions for the Poisson solve, including the zonal average;
  • option checks;
  • verification: ion-acoustic Landau damping along B0, quasi-neutral (gamma -0.168 vs -0.153) and with the vacuum term (-0.179 vs -0.194).
• Example slab_drift_wave: a drift wave in a density gradient.
  • Equilibrium: the initial state must be the exact Vlasov equilibrium f0(x + eps v_y), here a Maxwellian shifted by the ion diamagnetic velocity.
  • Time step: in the quasi-neutral limit the explicit splitting needs dt c_s pi/h <~ 2, since the ion-acoustic branch has no Debye cutoff.
  • Result: measured omega = 0.262, against 0.262 from the full-orbit dispersion relation and 0.287 from the drift-kinetic one (the O(k_perp^2 rho_i^2) gap).

Example slab_drift_wave for DriftKineticElectrostaticAdiabatic

The same physical setup as the full-orbit example (box, eps = 0.1, kappa_n = 1, T_i = T_e = 1, mode k = (pi/2, 3, 0.05)), for a direct comparison:

• drift-kinetic ions with the guiding-center Maxwellian n0(X) M (an exact equilibrium);
• explicit guiding-center pushers, dt = 0.25;
• LoadingParameters and GyroMaxwellian2D with B0 = 1 (both default to 2);
• Same electron model as the full-orbit example: Boltzmann response for all modes, i.e. no flux-surface average. For this, GyrokineticPoissonSolve gets the new option which_geometry=None; the default stays "cylindrical".
  • In the shearless slab the physical choice would be the average along the field lines (over z, "cylindrical"), which removes the electron response of the k_z = 0 modes.
  • That is not possible in the quasi-neutral full-orbit model, where the field equation has no polarization term to determine , and it is not needed for the studied mode, which has k_z != 0.
  • With "cylindrical", the ion convective cells of the k_z = 0 modes (omega ~ 3) also made the explicit push unstable at dt = 0.25.
• Result: measured omega = 0.291 against 0.287 from the drift-kinetic dispersion relation; the full-orbit example gives 0.262 for a predicted 0.262. pproc_slab_drift_wave.py --compare <full-orbit run> overlays both modes.

Fixes in the full-f Poisson right-hand side (separate commits)

1. ImplicitDiffusion: the full-f neutralizing background is -n0, not -Z n0.
   • The H1 kernels deposit the particle density, and the charge (number and sign) is in the coefficient (alpha^2/eps, with eps ~ 1/Z).
   • With -Z n0, the right-hand side was n - Z n0, which is wrong for Z != 1. Electrons (e.g. ColdPlasmaVlasov) got n + n0.
   • All existing examples and tests use Z = 1.
2. Initial Poisson solve: delta-f deposit without double subtraction.
   • VlasovAmpereOneSpecies, VlasovMaxwellOneSpecies and ColdPlasmaVlasov deposited the delta-f weights for the initial solve, and ImplicitDiffusion added the full-f background on top: n - 2 n0. In periodic examples the compatibility projection hid the offset; with Dirichlet BCs or a non-uniform n0 the initial phi was wrong.
   • New constructor argument ImplicitDiffusion(..., delta_f_source=True) (also PoissonSolve): particle sources are deposited with delta-f weights also for full-f markers (weights restored afterwards), and no background source is added. The three models use it for initial_poisson.
   • The delta-f deposit has far less noise than full-f weights minus n0: for a 1 % perturbation the initial phi has a relative error of 1 % instead of 400 % (100 markers per cell).
   • DeltaFParticles6D (LinearVlasovAmpere/MaxwellOneSpecies) no longer gets the -n0 background: their weights are delta-f already.
   • VlasovMaxwellOneSpecies gauss_error for full-f markers: the diagnostic deposits the full-f weights (whose charge the time stepping conserves) minus a cached t = 0 offset alpha^2/eps (rho_full(0) - rho_delta_f(0)). After 20 steps the residual is about 0.1 of the field term instead of about 3 times it (f0 deposit noise).
3. GyrokineticPoissonSolve docstrings corrected:
   • stab_mat default "M0ad_withT", not "Id";
   • the solved equation includes phi - <phi>, 1/eps and the rho_coeffs;
   • the parameters are those of the constructor.

Each of fixes 1 and 2 comes with tests that fail on the old code: phi_rms 2.0 (Z = 2 background), 0.09 (double subtraction, background for delta-f markers), relative error 4.1 (full-f deposit).

Tests (local)

• test_adiabatic_poisson.py (serial, 2 ranks); test_poisson.py, test_gyrokinetic_poisson.py (2 ranks, 140 passed)
• verification tests of VlasovPoissonElectrostaticAdiabatic, VlasovAmpereOneSpecies, VlasovMaxwellOneSpecies (2 ranks)
• default-params model tests of the new model, DriftKineticElectrostaticAdiabatic, VlasovAmpereOneSpecies, VlasovMaxwellOneSpecies, ColdPlasmaVlasov (4 ranks)
• example --test-run + pproc, ruff, doc/check_models_doc.py, doc/check_propagator_doc.py

Follow-ups (not in this PR)

• In the slab drift-wave example the kinetic energy of f - f0 drifts by 1-2 % of the thermal energy over t = 15 (unexplained).

Related: #799 (propagator order of DriftKineticElectrostaticAdiabatic, merged).

🤖 Generated with Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

An organization owner can choose the repository organization as the billing source for reviews requested by members with a paid Copilot license. Learn about Copilot code review billing.