Fix total energy at t=0 in Vlasov-Ampère/Maxwell models (ghost regions, control variate)

[max-models]

Problem

Found while writing the GPU end-to-end test of VlasovAmpereOneSpecies (#734). Setup: VlasovAmpereOneSpecies(alpha=1, epsilon=-1, with_B0=False), Cuboid r1=12.56, 8x1x1 elements, degree (3,1,1), ppc=40 sobol_standard, Maxwellian background, ModesCos density perturbation of amplitude 0.1, dt=0.1, 3 steps, solver tolerances 1e-12, initial Poisson with stab_eps=1e-6.

total_energy (electric_energy + kinetic_energy):

	t=0	step 1	step 2	step 3
full-f, before	18.96461097	18.96766028	18.96766028	18.96766028
full-f, after	18.96461097	18.96461097	18.96461097	18.96461097
control variate, before	18.96461097	0.11536371	0.09705431	0.08242210
control variate, after	0.12553637	0.10389888	0.08562030	0.07105039

Full-f: the Crank-Nicolson coupling conserved the total to round-off from step 1 on, but step 1 itself jumped by 1.6e-4 relative. Control variate: the t=0 value was a different quantity than the later values.

Root causes and fix

1. Ghost regions of the initial E field. post_allocate sets e = -grad(phi) with derham.grad.dot(-phi, out=e), which does not fill the ghost regions. The pusher in VlasovAmpereCoupling evaluates e_mid = (e^n + e^{n+1})/2 at the markers, ghost cells included, so in the first step it pushed with a wrong field near the periodic boundary. The field solve itself does not read ghost cells, so the field and the velocities were no longer consistent and energy was not conserved. From step 2 on, e comes from update_feec_variables, which syncs the ghost regions. This changes the physics of the first step, not only the diagnostic. Fix: e.update_ghost_regions() after the initial Poisson solve in VlasovAmpereOneSpecies, VlasovMaxwellOneSpecies, ColdPlasmaVlasov and LinearVlasovAmpereOneSpecies. These are all the models that write -grad(phi) into e_field this way.

2. Control-variate weights reset after the initial Poisson solve. initialize_weights already sets delta-f weights when control_variate=True. post_allocate calls update_weights() so that the Poisson right-hand side is the charge of f - f0, and then set weights = weights0 without checking the control variate. That left full-f weights for the t=0 diagnostic. The first PushEta then called update_weights(), so from step 1 on the weights were delta-f again. Fix: reset to weights0 only without the control variate. The same change is applied to the VlasovMaxwellOneSpecies Gauss-law diagnostic, which restored full-f weights at every scalar update. ColdPlasmaVlasov had the opposite problem: it never reset, so its full-f runs used delta-f weights for the whole run. It now follows the same rule.

3. Sobol-loaded markers were not sent to their MPI process (pic/base.py). With loading="sobol_standard"/"sobol_antithetic" every process draws its markers on the whole unit cube, but mpi_sort_markers was only called when box sorting was enabled. Without it, the first accumulation (the initial Poisson solve) wrote outside the local stencil arrays and corrupted the heap. The new test_energy_conservation is the only MPI test that loads Sobol markers without box sorting, and it made the next test abort the 2-rank verification job on Linux (Fatal Python error: Aborted). On 2 ranks, half of each rank's 3200 markers were outside its domain before the fix and none after. Fix: call mpi_sort_markers() after a Sobol draw when there is an MPI communicator and box sorting is off.

With the control variate, kinetic_energy is now the kinetic energy of f - f0 at all times, including t=0. This is now stated in the KineticEnergyPIC docstring, and it matches what LinearMHDVlasovCC and ToyDrift already did.

Expected behaviour (no code change)

With the control variate, the delta-f total energy is not conserved to round-off, even after this fix. The Crank-Nicolson step conserves E_E + α²/2 Σ w_δ |v|² with the weights frozen during the step. Afterwards update_weights changes w_δ by -Δf0(v)/(s0 N), which adds two things per step:

• a Monte-Carlo estimate of zero (∝ Σ f0 (v·E)|v|²), and
• an O(Δt²) term: the energy the field gives to the background markers.

Measured on 16 elements with ppc=4000 up to T=0.3, the drift is 3.0e-3, 1.2e-3 and 2.7e-4 for dt = 0.1, 0.05 and 0.025. So it converges with dt, and it is a property of the scheme, not of the diagnostic.

Affected models

VlasovAmpereOneSpecies, VlasovMaxwellOneSpecies and ColdPlasmaVlasov (ghost regions + weights), LinearVlasovAmpereOneSpecies (ghost regions only; its markers are delta-f by construction).

Tests

• New test_energy_conservation[False/True] in models/tests/verification/test_verif_VlasovAmpereOneSpecies.py. It uses the setup above with 16 elements, so that the CI run on 4 ranks with clones gets at least p=3 elements per rank. ~2 s each.
  • Full-f, ppc=40: drift from t=0 must be < 1e-10. On devel it is 9.4e-7; after the fix it is round-off.
  • Control variate, ppc=400: drift must be < 0.1, for the reason above. On devel it is 0.99; after the fix it is 3.7e-3.
  • Both fail on devel and pass with the fix: serial, mpirun -n 2 --with-mpi, mpirun -n 4 --with-mpi and mpirun -n 4 --with-mpi --nclones 2.
• Passing: test_weak_Landau (control variate), test_verif_LinearVlasovAmpereOneSpecies.py, test_verif_VlasovMaxwellOneSpecies.py, and the default-parameter runs (test_single_model) of VlasovAmpereOneSpecies, VlasovMaxwellOneSpecies and LinearVlasovAmpereOneSpecies.
• Not run to completion locally: the default-parameter ColdPlasmaVlasov run takes more than 5 minutes on devel as well. I stopped it, so CI covers it.

Note for #734

The GPU e2e test of VlasovAmpereOneSpecies in #734 currently uses full-f weights and only a loose t=0 check (|en_tot[1] - en_tot[0]| < 1e-3) because of these two bugs. Once this is merged, it can check total_energy from t=0 at solver tolerance, and it can also use the control variate.

🤖 Generated with Claude Code