Initial Poisson solve: default stabilization puts net charge / 1e-14 into phi, E loses digits #766

Bug: the initial Poisson solve of VlasovAmpereOneSpecies (and other models that use PoissonSolve with the defaults) loses digits in E = -grad(phi), because of how the periodic problem is stabilized.

PoissonSolve.Options.stab_eps defaults to 0.0, and ImplicitDiffusion then raises it to sigma_1 = 1e-14 (propagators/implicit_diffusion.py:289). With periodic boundary conditions, the Monte-Carlo charge density of a finite number of markers has a small nonzero net charge. The stabilized solve puts (net charge) / 1e-14 into the constant mode of phi. So phi is about 1e14, and E = -grad(phi) loses about 14 digits to cancellation.

Observed (reported in #734, Landau damping setup with 320 Sobol markers): about 5% error in E, on NumPy and on CuPy. test_verif_VlasovAmpereOneSpecies.py uses this default (Options(stab_mat="M0"), line 81). The LinearVlasovAmpere verification tests and the e2e GPU tests of #734 set stab_eps=1e-6, which hides the issue.

Possible fixes (to discuss):

• remove the constant mode from the right-hand side (or project it out of phi) for periodic problems, since only grad(phi) is used;
• or a larger default stabilization, documented together with its effect on the constant mode;
• in any case, use the same setting in the verification tests of both VlasovAmpere models.

🤖 Generated with Claude Code