Repository navigation
Rename VlasovPoissonElectrostaticAdiabatic to VlasovPoissonAdiabatic - #804
Merged
Merged
Conversation
spossann
added this pull request to stack #806
October 20, 2026 14:34
The Poisson equation already says that the model is electrostatic. Renames the model module, the class, the example folder and the verification test, and updates the references in the docs, the API and the drift-kinetic slab example. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
spossann
force-pushed
the
rename-vp-adiabatic
branch
from
October 20, 2026 15:17
d857565 to
922f3da
Compare
spossann
marked this pull request as ready for review
October 20, 2026 15:17
spossann
added a commit
that referenced
this pull request
Oct 21, 2026
…ssical Pauli particles (#805) New model `VlasovPoissonAdiabaticPauli`: guiding-center ions computed with **classical Pauli particles** (J. Xiao and H. Qin, "Slow manifolds of classical Pauli particle enable structure-preserving geometric algorithms for guiding center dynamics", Comput. Phys. Commun. 265, 107981 (2021)). It uses the existing full-orbit pushers, with time steps larger than the cyclotron period. Stacked on #804 (rename of `VlasovPoissonElectrostaticAdiabatic` to `VlasovPoissonAdiabatic`). Merge that first, then retarget this PR to `devel`. ## What's in it - **`PauliParticles6D`** (`pic/particles.py`): full-orbit markers with a magnetic moment in one diagnostics column. - After the usual drawing and weight initialization, each marker gets mu = |v_perp|^2/(2|B0|), and its velocity is set to v -> (v.b0) b0. This is Eq. (7) of the paper and the paragraph after it: the marker starts on the slow manifold, at its guiding-center position. - With the control variate, f0 is evaluated at the gyrotropic velocity v_par b0 + sqrt(2 mu |B0|) e_perp. The backgrounds must therefore be isotropic Maxwellians. - `mu_prime()` returns the paper's mu' = |v x b0|^2/(2|B0|). - Markers loaded from restart or external files keep their mu. - **`PushVinForceField(magnetic_moment=True)`** with the new kernel `push_v_with_efield_mu`. - The kick uses E† = -grad(phi + epsilon mu |B0|) (Eq. (12)), with the same projected grad|B0| (`gradB1`) as the guiding-center pushers. - `PushVxB` is unchanged. With Strang splitting and `algo="implicit"` (Cayley rotation), one step is the Boris scheme of the paper (BAP2). - **`AdiabaticPoissonSolve.Options(polarization=True)`** (quasi-neutral only): adds the guiding-center polarization term -div(n0/|B0|^2 grad_perp phi) (`M1gyro`), the same operator as `GyrokineticPoissonSolve`. - **`VlasovPoissonAdiabaticPauli`** (subclass of `VlasovPoissonAdiabatic`): - `with_poisson=False` for orbit studies in the equilibrium field; - scalars `en_mu`, `mu_prime_max`, and the polarization energy in `en_phi`; - a warning for delta-f without the polarization term, and for full-f with it. - **Examples**: `tokamak_orbits` (banana and passing orbits, dt = 8 cyclotron periods) and `slab_drift_wave` (delta-f with the polarization term, the same setup as the other two Poisson models). - Fix of the CI failure: the kernel count in `test_catalog_signatures` is now 45. ## Polarization with delta-f The marker positions are guiding-center positions, and on the slow manifold the polarization drift compresses the markers. - **Full-f:** this compression already gives the polarization density. - **Delta-f:** the control variate subtracts f0 along the markers, and this removes the compression from the deposited charge. So the polarization term belongs in the quasi-neutrality equation, as for guiding centers. Slab drift wave (delta-f, 8x8x4 elements, 100 ppc, 2 ranks): | dt / cyclotron period | polarization term | omega (DK / Pauli prediction 0.287) | max mu' at t = 25 | |---|---|---|---| | 0.08 | no | 0.316 (DK root without polarization: 0.317) | blows up from t ~ 10 | | 8 | no | 0.315 | blows up from t ~ 10 | | 0.08 | yes | 0.291 | 4.4e-3 | | 8 | yes | 0.292 | 3.6e-3 | | 16 | yes | 0.294 | 4.6e-3 | **Open issue.** A weak instability remains with the polarization term. It needs the density gradient and is independent of dt: mu' grows by about x5 per 5 time units, and the field energy grows from t ~ 25 (about one wave period) to saturation at t ~ 40. `DriftKineticElectrostaticAdiabatic` stays quiet to t = 50 in the same setup. - Removing the advection part n0'·grad phi from the polarization operator, which the delta-f weights also carry, does not change it (tested with a non-symmetric operator). - A likely cause is that the compression part is now instantaneous (field equation), while the advection part follows the marker orbits, including their cyclotron response. - The example runs to t = 25, where the drift wave is clean. The model documentation describes the issue. ## Orbit verification `test_verif_VlasovPoissonAdiabaticPauli.py` compares Pauli orbits with the explicit RK4 pushers of `GuidingCenter` (dt = 0.8 cyclotron periods) in an ad hoc tokamak, epsilon = 0.01, for trapped and passing markers: | dt / cyclotron period | distance from GC orbit (passing) | energy error | max mu'/mu | |---|---|---|---| | 16 | 0.050 - 0.064 | 6.2e-3 | 6.8e-4 | | 8 | 0.012 - 0.016 | 1.6e-3 | 2.0e-4 | | 4 | 0.003 - 0.004 | 4.1e-4 | 2.3e-4 | | 0.8 | 0.0005 | 2.1e-5 | 3.0e-4 | The Pauli orbits converge to the GC orbits at second order in dt (orbit widths are about 0.04 - 0.09). ## Tests - `pic/tests/test_pauli_particles.py`: the same markers, sampling densities and weights (full-f and control variate) as `Particles6D`, in a tilted uniform field and in a screw pinch. - `propagators/tests/test_adiabatic_poisson.py`: `AdiabaticPoissonSolve(polarization=True)` against the analytic solution (as `GyrokineticPoissonSolve`), plus option checks. - The orbit verification test, serially and on 4 ranks; the default-params model tests on 4 ranks; the docs checks; both examples with `--test-run`. Side note (not changed here): a fancy-indexed positions array of a single marker has strides (8, 24) and is still flagged C-contiguous. So `ascontiguousarray` in the marker evaluation of `Domain` keeps it, and the Pyccel kernels reject it. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
spossann
added a commit
that referenced
this pull request
Oct 21, 2026
…804) Renames the model `VlasovPoissonElectrostaticAdiabatic` (added in #800) to `VlasovPoissonAdiabatic`: the Poisson equation already says that the model is electrostatic. - model module `vlasov_poisson_adiabatic.py`, class, example folder `examples/VlasovPoissonAdiabatic/` and verification test `test_verif_VlasovPoissonAdiabatic.py` - references in `models-all.rst`, `struphy.api.models`, the generated `models/__init__.py` and the drift-kinetic slab example - no deprecated alias (the old name was never released) Tests: the verification test, the default-params model test on 4 ranks, `doc/check_models_doc.py` and both slab drift-wave examples (and pprocs) with `--test-run`. The new model `VlasovPoissonAdiabaticPauli` is stacked on this branch: #805. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
spossann
added a commit
that referenced
this pull request
Oct 21, 2026
…ssical Pauli particles (#805) New model `VlasovPoissonAdiabaticPauli`: guiding-center ions computed with **classical Pauli particles** (J. Xiao and H. Qin, "Slow manifolds of classical Pauli particle enable structure-preserving geometric algorithms for guiding center dynamics", Comput. Phys. Commun. 265, 107981 (2021)). It uses the existing full-orbit pushers, with time steps larger than the cyclotron period. Stacked on #804 (rename of `VlasovPoissonElectrostaticAdiabatic` to `VlasovPoissonAdiabatic`). Merge that first, then retarget this PR to `devel`. ## What's in it - **`PauliParticles6D`** (`pic/particles.py`): full-orbit markers with a magnetic moment in one diagnostics column. - After the usual drawing and weight initialization, each marker gets mu = |v_perp|^2/(2|B0|), and its velocity is set to v -> (v.b0) b0. This is Eq. (7) of the paper and the paragraph after it: the marker starts on the slow manifold, at its guiding-center position. - With the control variate, f0 is evaluated at the gyrotropic velocity v_par b0 + sqrt(2 mu |B0|) e_perp. The backgrounds must therefore be isotropic Maxwellians. - `mu_prime()` returns the paper's mu' = |v x b0|^2/(2|B0|). - Markers loaded from restart or external files keep their mu. - **`PushVinForceField(magnetic_moment=True)`** with the new kernel `push_v_with_efield_mu`. - The kick uses E† = -grad(phi + epsilon mu |B0|) (Eq. (12)), with the same projected grad|B0| (`gradB1`) as the guiding-center pushers. - `PushVxB` is unchanged. With Strang splitting and `algo="implicit"` (Cayley rotation), one step is the Boris scheme of the paper (BAP2). - **`AdiabaticPoissonSolve.Options(polarization=True)`** (quasi-neutral only): adds the guiding-center polarization term -div(n0/|B0|^2 grad_perp phi) (`M1gyro`), the same operator as `GyrokineticPoissonSolve`. - **`VlasovPoissonAdiabaticPauli`** (subclass of `VlasovPoissonAdiabatic`): - `with_poisson=False` for orbit studies in the equilibrium field; - scalars `en_mu`, `mu_prime_max`, and the polarization energy in `en_phi`; - a warning for delta-f without the polarization term, and for full-f with it. - **Examples**: `tokamak_orbits` (banana and passing orbits, dt = 8 cyclotron periods) and `slab_drift_wave` (delta-f with the polarization term, the same setup as the other two Poisson models). - Fix of the CI failure: the kernel count in `test_catalog_signatures` is now 45. ## Polarization with delta-f The marker positions are guiding-center positions, and on the slow manifold the polarization drift compresses the markers. - **Full-f:** this compression already gives the polarization density. - **Delta-f:** the control variate subtracts f0 along the markers, and this removes the compression from the deposited charge. So the polarization term belongs in the quasi-neutrality equation, as for guiding centers. Slab drift wave (delta-f, 8x8x4 elements, 100 ppc, 2 ranks): | dt / cyclotron period | polarization term | omega (DK / Pauli prediction 0.287) | max mu' at t = 25 | |---|---|---|---| | 0.08 | no | 0.316 (DK root without polarization: 0.317) | blows up from t ~ 10 | | 8 | no | 0.315 | blows up from t ~ 10 | | 0.08 | yes | 0.291 | 4.4e-3 | | 8 | yes | 0.292 | 3.6e-3 | | 16 | yes | 0.294 | 4.6e-3 | **Open issue.** A weak instability remains with the polarization term. It needs the density gradient and is independent of dt: mu' grows by about x5 per 5 time units, and the field energy grows from t ~ 25 (about one wave period) to saturation at t ~ 40. `DriftKineticElectrostaticAdiabatic` stays quiet to t = 50 in the same setup. - Removing the advection part n0'·grad phi from the polarization operator, which the delta-f weights also carry, does not change it (tested with a non-symmetric operator). - A likely cause is that the compression part is now instantaneous (field equation), while the advection part follows the marker orbits, including their cyclotron response. - The example runs to t = 25, where the drift wave is clean. The model documentation describes the issue. ## Orbit verification `test_verif_VlasovPoissonAdiabaticPauli.py` compares Pauli orbits with the explicit RK4 pushers of `GuidingCenter` (dt = 0.8 cyclotron periods) in an ad hoc tokamak, epsilon = 0.01, for trapped and passing markers: | dt / cyclotron period | distance from GC orbit (passing) | energy error | max mu'/mu | |---|---|---|---| | 16 | 0.050 - 0.064 | 6.2e-3 | 6.8e-4 | | 8 | 0.012 - 0.016 | 1.6e-3 | 2.0e-4 | | 4 | 0.003 - 0.004 | 4.1e-4 | 2.3e-4 | | 0.8 | 0.0005 | 2.1e-5 | 3.0e-4 | The Pauli orbits converge to the GC orbits at second order in dt (orbit widths are about 0.04 - 0.09). ## Tests - `pic/tests/test_pauli_particles.py`: the same markers, sampling densities and weights (full-f and control variate) as `Particles6D`, in a tilted uniform field and in a screw pinch. - `propagators/tests/test_adiabatic_poisson.py`: `AdiabaticPoissonSolve(polarization=True)` against the analytic solution (as `GyrokineticPoissonSolve`), plus option checks. - The orbit verification test, serially and on 4 ranks; the default-params model tests on 4 ranks; the docs checks; both examples with `--test-run`. Side note (not changed here): a fancy-indexed positions array of a single marker has strides (8, 24) and is still flagged C-contiguous. So `ascontiguousarray` in the marker evaluation of `Domain` keeps it, and the Pyccel kernels reject it. 🤖 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
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.
Renames the model
VlasovPoissonElectrostaticAdiabatic(added in #800) toVlasovPoissonAdiabatic: the Poisson equation already says that the model is electrostatic.vlasov_poisson_adiabatic.py, class, example folderexamples/VlasovPoissonAdiabatic/and verification testtest_verif_VlasovPoissonAdiabatic.pymodels-all.rst,struphy.api.models, the generatedmodels/__init__.pyand the drift-kinetic slab exampleTests: the verification test, the default-params model test on 4 ranks,
doc/check_models_doc.pyand both slab drift-wave examples (and pprocs) with--test-run.The new model
VlasovPoissonAdiabaticPauliis stacked on this branch: #805.🤖 Generated with Claude Code