Repository navigation
DriftKineticElectrostaticAdiabatic: Poisson solve between the two guiding-center pushers - #799
Merged
Merged
Conversation
…ding-center pushers Reorder the propagators from [gc_poisson, push_gc_bxe, push_gc_para] to [push_gc_bxe, gc_poisson, push_gc_para]. Both pushers move the guiding centers, so with Strang splitting the old order pushed para(dt) and the closing bxe(dt/2) with a potential solved before a drift. Now the Strang step is bxe(dt/2), P, para(dt), P, bxe(dt/2) at the same cost of two Poisson solves. Also solve the Poisson equation in update_scalar_quantities (once per saved step), so that en_phi, en_tot and the saved phi belong to the current marker positions. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
spossann
commented
Oct 19, 2026
… data Remove the update_scalar_quantities override from the previous commit (review of #799): solving the Poisson equation once per saved step is too expensive. The saved potential and en_phi belong to the marker positions before the last push of a time step; this is stated in doc_discretization. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
spossann
marked this pull request as ready for review
October 20, 2026 06:22
spossann
added a commit
that referenced
this pull request
Oct 20, 2026
…atic electrons) and fixes of the full-f Poisson background (#800) ## 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 #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 <phi>, 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](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Both guiding-center pushers of
DriftKineticElectrostaticAdiabaticmove the markers and use the potential. With the old propagator order[gc_poisson, push_gc_bxe, push_gc_para], a Strang step (StruphyModel.integrate: all but the last with dt/2, the last with dt, then the others reversed) wasso
para(dt)and the closingbxe(dt/2)pushed with a potential solved before a drift; the two back-to-back Poisson solves at the step boundary were redundant.This PR reorders to
[push_gc_bxe, gc_poisson, push_gc_para]:para(dt)and the closingbxe(dt/2)use a freshly solved potential; the twobxe(dt/2)at a step boundary act as onebxe(dt)with the potential solved before it.phianden_phibelong to the marker positions before that last push (bxe(dt/2)with Strang,para(dt)with Lie-Trotter); this is stated indoc_discretization.doc_discretizationexplains the order; the class link now points toGyrokineticPoissonSolve.Note on Lie-Trotter
Both examples (
itg_cylindre,cyclone) usesplit_algo="LieTrotter". There the reorder only moves the stale potential fromparatobxe(one push per step uses the potential of the previous step in either order), so their results may change slightly without becoming more accurate. The gain is for Strang splitting.Tests
pytest -m single --model-name DriftKineticElectrostaticAdiabatic src/struphy/models/tests/default_params/(serial)--test-rundoc/check_models_doc.pyFound while implementing the full-orbit counterpart
VlasovPoissonElectrostaticAdiabatic(separate branchvlasov-poisson-adiabatic, not yet pushed), which uses the order drift, Poisson, E-kick, v x B rotation for the same reason.🤖 Generated with Claude Code