Remove net charge from the rhs of singular (periodic) Poisson solves

[max-models]

Problem

On a domain without Dirichlet boundary conditions (e.g. fully periodic), the Poisson operator G^T M1 G is singular. Its kernel is the constant functions, which have coefficient vector 1 (all ones, by partition of unity). PoissonSolve only regularized this with stab_eps * stab_mat, and the default stab_eps = 0 is bumped to 1e-14 in allocate. When the source has a nonzero net charge 1^T b, the constant mode of phi becomes (net charge) / stab_eps. In PIC this always happens, because the Monte-Carlo charge density is noisy. That constant is around 1e9 to 1e13, so E = -grad phi loses most of its digits to cancellation.

This affects the initial_poisson of VlasovAmpereOneSpecies, VlasovMaxwellOneSpecies, LinearVlasovAmpereOneSpecies, LinearVlasovMaxwellOneSpecies and ColdPlasmaVlasov. It also affects every other PoissonSolve user on periodic or Neumann domains (Poisson, HasegawaWakatani, IncompressibleNavierStokesSPH, ToyDrift). The verification tests use the default stab_eps, so they are affected too.

I found this while writing the GPU end-to-end test for VlasovAmpereOneSpecies (#734). That test currently sets stab_eps=1e-6 as a workaround, which can be removed once this PR is merged.

Numbers before the fix

All runs use 1D periodic Cuboid(r1=12.56), degree (3,1,1) and default options.

setup	max|phi|	rel. error of E
analytic source 0.3 + 0.1 cos(kx), 16 cells, stab_mat="Id"	2.4e13	4.6e-2
same, stab_mat="M0"	3.0e13	2.7e-2
same, stab_eps=1e-6 (workaround)	3.0e5	4.8e-4
PIC, control variate, pseudo_random, ppc 40, 8 cells (net charge 7.9e-2)	6.3e11	differs by 2.7e-4 from the stab_eps=1e-6 result
test_verif_VlasovAmpereOneSpecies initial solve	3.8e8 (coefficients)

After the fix

setup	max|phi|	rel. error of E
analytic source, Id or M0, default stab_eps	0.40 (exact: 0.40)	4.8e-4 (discretization error)
PIC case above	0.48	differs by 3.9e-6 from stab_eps=1e-6 (that difference is the O(stab_eps) error of the workaround itself)
test_verif_VlasovAmpereOneSpecies initial solve	4.0e-3 (coefficients)

Fix

New option enforce_compatibility. It is True for PoissonSolve and False for ImplicitDiffusion and PoissonAdiabaticGyrokinetic, where the mass term is physical. When it is on and the operator is singular (no Dirichlet BCs and no polar splines, see ImplicitDiffusion.operator_is_singular), the rhs is made to satisfy the discrete compatibility condition 1^T b = 0 before the solve:

b <- b - (1^T b) / (1^T M0 1) * M0 1

• Physics: M0 1 is the weak form of a uniform density, so this adds a uniform neutralizing background. A periodic plasma is quasi-neutral, and the net Monte-Carlo charge is noise that must be discarded.
• Effect on phi: phi stays O(1). stab_eps now only fixes the constant of phi, which is harmless, so its default is unchanged. With stab_mat="M0" the projection changes phi only by a constant for any stab_eps, so E is the same; the new test checks this.
• Implementation: it uses only feectools vector operations (inner, which is a global MPI reduction, and mul_iadd), with 1 and M0 1 precomputed in allocate. There are no host copies, so it should work on both cunumpy backends. The CuPy path is untested.
• Out of scope: polar splines are conservatively treated as regular and left unchanged, because the polar coefficient vector of the constant function is not simply all ones.
• Side effect: CG converges faster on the singular problems. In test_poisson_2d_multigrid (periodic, unpreconditioned) it now takes 97 iterations instead of 142.

Tests

• New test_poisson_net_charge_periodic_1d[Id|M0] (about 1 s). It uses a non-neutral analytic source with default options and checks:
  • phi is O(1) and matches the exact solution up to a constant;
  • the E error is below 1e-3;
  • with enforce_compatibility=False, |phi| > 1e10;
  • for M0 with finite stab_eps, the projection only shifts phi by a constant.
• New test_poisson_dirichlet_not_singular: with Dirichlet BCs the rhs is left unchanged.
• The new tests also pass with mpirun -n 2.
• Locally passing: the 40 serial tests of test_poisson.py before the multigrid ones; test_poisson_2d_multigrid[periodic-degree0] (about 3 min per run, so the rest of the multigrid set is left to CI); test_implicit_diffusion_multigrid_dt; the 24 test_poisson_M1perp_1d cases in test_gyrokinetic_poisson.py; test_verif_VlasovAmpereOneSpecies.py (48 s).
• The 2D/3D gyrokinetic Poisson tests and the remaining verification tests are left to CI.

🤖 Generated with Claude Code