From 8e5ed8a77f5c5be1a1e42d1d8be77b442a22c5b9 Mon Sep 17 00:00:00 2001 From: Max Date: Thu, 8 Oct 2026 23:00:42 +0200 Subject: [PATCH] Remove net charge from rhs of singular Poisson solves On a periodic domain (no Dirichlet BCs) the Poisson operator G^T M1 G is singular with the constant functions (coefficient vector of all ones) in its kernel. It was regularized only with stab_eps (default -> 1e-14), so a non-neutral source (e.g. Monte-Carlo noise of the particle charge density) made the constant mode of phi ~ (net charge) / stab_eps ~ 1e13, and E = -grad phi lost most digits to cancellation (2-5 % error). PoissonSolve now enforces the discrete compatibility condition 1^T b = 0 by removing the net charge from the rhs, b <- b - (1^T b)/(1^T M0 1) M0 1 (a uniform neutralizing background), when the operator is singular (no Dirichlet BCs, no polar splines). New option enforce_compatibility (True for PoissonSolve, False for ImplicitDiffusion and PoissonAdiabaticGyrokinetic). Co-Authored-By: Claude Opus 5.5 --- src/struphy/propagators/implicit_diffusion.py | 45 ++++++++- .../poisson_adiabatic_gyrokinetic.py | 2 + src/struphy/propagators/poisson_solve.py | 20 +++- src/struphy/propagators/tests/test_poisson.py | 93 +++++++++++++++++++ 4 files changed, 158 insertions(+), 2 deletions(-) diff --git a/src/struphy/propagators/implicit_diffusion.py b/src/struphy/propagators/implicit_diffusion.py index 2af8fcd03..1da467b52 100644 --- a/src/struphy/propagators/implicit_diffusion.py +++ b/src/struphy/propagators/implicit_diffusion.py @@ -134,7 +134,7 @@ def __init__( @dataclass(repr=False) class Options(OptionsBase): - """Configuration options for :class:`ImplicitDiffusion`. + r"""Configuration options for :class:`ImplicitDiffusion`. Parameters ---------- @@ -199,6 +199,22 @@ class Options(OptionsBase): ``verbose``, ``info``, ``recycle``). If ``None``, defaults to ``SolverParameters()``. + enforce_compatibility : bool, default=False + If ``True`` and the diffusion operator is singular (no Dirichlet boundary conditions + and no polar splines, i.e. the constant functions are in its kernel), the right-hand side + is made compatible before the solve by removing its constant ("net charge") component, + + .. math:: + + \mathbf b \leftarrow \mathbf b - \frac{\mathbf 1^\top \mathbf b}{\mathbf 1^\top \mathbb M^0 \mathbf 1}\, \mathbb M^0 \mathbf 1\,, + + where :math:`\mathbf 1` (all ones) is the coefficient vector of the constant function + (partition of unity). This is equivalent to adding a uniform neutralizing background. + It is meant for (nearly) singular problems where ``sigma_1`` is only a small + stabilization; then ``phi`` stays of order one instead of growing like + ``(net charge) / sigma_1``, which would destroy the accuracy of ``grad phi``. + For ``stab_mat="M0"`` and ``sigma_2 = 0`` it changes ``phi`` only by a constant, for any ``sigma_1``. + filter_params : dict[PICVariable | SPHVariable, FilterParameters], default=None If not None, specifies a filter to the accumulation of a specific variable. Keyed by the ``pic_variable`` of the corresponding @@ -226,6 +242,7 @@ class Options(OptionsBase): precond: LiteralOptions.OptsDiffusionPrecond = "MassMatrixPreconditioner" multigrid: MultiGridOptions = None solver_params: SolverParameters = None + enforce_compatibility: bool = False filter_params: dict[PICVariable | SPHVariable, FilterParameters] = None def __post_init__(self): @@ -391,6 +408,28 @@ def verify_rhs(rho) -> StencilVector | FEECVariable | AccumulatorVector: self._rhs2 = phi.space.zeros() self._tmp_src = phi.space.zeros() + # constant function and its weak form M0 @ 1, for removing the kernel component of the rhs + self._ones = None + if self.options.enforce_compatibility and self.operator_is_singular: + self._ones = phi.space.zeros() + owned = tuple(slice(s, e + 1) for s, e in zip(self._ones.starts, self._ones.ends)) + self._ones[owned] = 1.0 + self._ones.update_ghost_regions() + self._m0_ones = self.mass_ops.M0.dot(self._ones) + self._ones_m0_ones = self._ones.inner(self._m0_ones) + + @property + def operator_is_singular(self) -> bool: + """True if the diffusion operator ``grad.T @ diffusion_mat @ grad`` is known to be singular + with the constant functions (coefficient vector of all ones) in its kernel: no Dirichlet boundary + conditions and tensor-product splines. Polar splines are conservatively treated as regular.""" + no_dirichlet = not any(any(bc) for bc in self.derham.dirichlet_bc) + return ( + no_dirichlet + and not self.derham.polar_splines + and isinstance(self.variables.phi.spline.vector, StencilVector) + ) + @property def sources(self) -> list[StencilVector | FEECVariable | AccumulatorVector]: """ @@ -475,6 +514,10 @@ def __call__(self, dt): rhs += self._rhs2 + # remove the kernel component (net charge) of the rhs: 1^T rhs = 0 (compatibility condition) + if self._ones is not None: + rhs.mul_iadd(-self._ones.inner(rhs) / self._ones_m0_ones, self._m0_ones) + if self.diagnostic is not None: proj = L2Projector("H1", self.mass_ops) self.diagnostic.spline.vector = proj.solve(rhs) diff --git a/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py b/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py index 7cacf5434..c1c5ea17b 100644 --- a/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py +++ b/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py @@ -165,6 +165,8 @@ class Options: precond: LiteralOptions.OptsMassPrecond = "MassMatrixPreconditioner" solver_params: SolverParameters = None filter_params: dict[PICVariable | SPHVariable, FilterParameters] = None + # the adiabatic term is physical, not a stabilization: the rhs is not projected + enforce_compatibility: bool = False def __post_init__(self): # checks diff --git a/src/struphy/propagators/poisson_solve.py b/src/struphy/propagators/poisson_solve.py index 9e71cb73c..59a0b9d93 100644 --- a/src/struphy/propagators/poisson_solve.py +++ b/src/struphy/propagators/poisson_solve.py @@ -35,7 +35,7 @@ class PoissonSolve(ImplicitDiffusion): @dataclass(repr=False) class Options(OptionsBase): - """Configuration options for :class:`Poisson`. + r"""Configuration options for :class:`Poisson`. Parameters ---------- @@ -81,6 +81,23 @@ class Options(OptionsBase): ``verbose``, ``info``, ``recycle``). If ``None``, defaults to ``SolverParameters()``. + enforce_compatibility : bool, default=True + If ``True`` and the Poisson operator is singular (no Dirichlet boundary conditions + and no polar splines, e.g. fully periodic), the constant ("net charge") component of the + right-hand side is removed before the solve, + + .. math:: + + \mathbf b \leftarrow \mathbf b - \frac{\mathbf 1^\top \mathbf b}{\mathbf 1^\top \mathbb M^0 \mathbf 1}\, \mathbb M^0 \mathbf 1\,, + + such that the discrete compatibility condition :math:`\mathbf 1^\top \mathbf b = 0` holds + (:math:`\mathbf 1` is the coefficient vector of the constant function). This is equivalent to + a uniform neutralizing background; for particle sources, the net charge is Monte-Carlo noise. + Without it, the constant mode of ``phi`` becomes ``(net charge) / stab_eps``, and + ``grad phi`` loses digits to cancellation. ``stab_eps`` then only fixes the constant of ``phi``. + Set to ``False`` if ``stab_eps * stab_mat`` is a physical term (e.g. adiabatic electrons with + ``stab_mat="M0ad"``) and not just a stabilization. + filter_params : dict[PICVariable | SPHVariable, FilterParameters], default=None If not None, specifies a filter to the accumulation of a specific variable. @@ -104,6 +121,7 @@ class Options(OptionsBase): precond: LiteralOptions.OptsDiffusionPrecond = "MassMatrixPreconditioner" multigrid: MultiGridOptions = None solver_params: SolverParameters = None + enforce_compatibility: bool = True filter_params: dict[PICVariable | SPHVariable, FilterParameters] = None def __post_init__(self): diff --git a/src/struphy/propagators/tests/test_poisson.py b/src/struphy/propagators/tests/test_poisson.py index 84524ee77..2484fbae7 100644 --- a/src/struphy/propagators/tests/test_poisson.py +++ b/src/struphy/propagators/tests/test_poisson.py @@ -880,6 +880,99 @@ def test_implicit_diffusion_multigrid_dt(): assert prop._solver._info["niter"] <= 15 +@pytest.mark.parametrize("stab_mat", ["Id", "M0"]) +def test_poisson_net_charge_periodic_1d(stab_mat): + """Periodic Poisson problem with a non-neutral source (net charge, e.g. Monte-Carlo noise). + + The constant functions span the kernel of the periodic Laplacian. With the default options the + net charge is removed from the rhs (compatibility condition), such that phi stays O(1) and + E = -grad phi is accurate. Without it, the constant mode of phi is (net charge) / stab_eps. + """ + Lx = 12.56 + k = 2 * xp.pi / Lx + amp = 0.1 + net = 0.3 + + domain = domains.Cuboid(r1=Lx) + derham = Derham(TensorProductGrid(num_elements=(16, 1, 1)), DerhamOptions(degree=(3, 1, 1)), comm=comm) + mass_ops = WeightedMassOperators(derham, domain) + Propagator.derham = derham + Propagator.domain = domain + Propagator.mass_ops = mass_ops + + def rho(e1, e2, e3): + return net + amp * xp.cos(2 * xp.pi * e1) + + # uniform points without the duplicate endpoint, such that the mean of cos(k x) vanishes + e1 = xp.linspace(0.0, 1.0, 100, endpoint=False) + x = domain(e1, 0.0, 0.0)[0][:, 0, 0] + phi_exact = amp / k**2 * xp.cos(k * x) + e_exact = amp / k * xp.sin(k * x) + + def solve(**kwargs): + phi = FEECVariable(space="H1") + phi.allocate(derham=derham, domain=domain) + poisson = PoissonSolve(rho=rho) + poisson.variables.phi = phi + poisson.options = poisson.Options( + stab_mat=stab_mat, + solver_params=SolverParameters(tol=1e-13, maxiter=3000, recycle=False), + **kwargs, + ) + poisson.allocate() + poisson(1.0) + e = FEECVariable(space="Hcurl") + e.allocate(derham=derham, domain=domain) + derham.grad.dot(-phi.spline.vector, out=e.spline.vector) + phi_vals = domain.push(phi.spline, e1, 0.0, 0.0, kind="0")[:, 0, 0] + e_vals = domain.push(e.spline, e1, 0.0, 0.0, kind="1")[0][:, 0, 0] + return poisson, phi_vals, e_vals + + # default options (stab_eps -> 1e-14): phi is O(1) and E is accurate + poisson, phi_vals, e_vals = solve() + assert poisson.operator_is_singular + e_err = xp.max(xp.abs(e_vals - e_exact)) / xp.max(xp.abs(e_exact)) + phi_err = xp.max(xp.abs(phi_vals - xp.mean(phi_vals) - phi_exact)) / xp.max(xp.abs(phi_exact)) + logger.info(f"{stab_mat=}: {xp.max(xp.abs(phi_vals))=}, {e_err=}, {phi_err=}") + assert xp.max(xp.abs(phi_vals)) < 2.0 * xp.max(xp.abs(phi_exact)) + assert e_err < 1e-3 + assert phi_err < 1e-3 + + # without compatibility, the constant mode of phi is huge: (net charge) / stab_eps + _, phi_vals_nc, _ = solve(enforce_compatibility=False) + logger.info(f"{stab_mat=}, enforce_compatibility=False: {xp.max(xp.abs(phi_vals_nc))=}") + assert xp.max(xp.abs(phi_vals_nc)) > 1e10 + + # for stab_mat="M0" and a finite stab_eps, the projection changes phi only by a constant + if stab_mat == "M0": + _, phi1, e1_vals = solve(stab_eps=1e-2) + _, phi2, e2_vals = solve(stab_eps=1e-2, enforce_compatibility=False) + assert xp.max(xp.abs(e1_vals - e2_vals)) < 1e-10 * xp.max(xp.abs(e_exact)) + shift = phi2 - phi1 + assert xp.max(xp.abs(shift - xp.mean(shift))) < 1e-10 * xp.max(xp.abs(phi_exact)) + + +def test_poisson_dirichlet_not_singular(): + """With Dirichlet boundary conditions the Poisson operator is regular and the rhs is not modified.""" + domain = domains.Cuboid(r1=4.0) + derham = Derham( + TensorProductGrid(num_elements=(8, 1, 1)), + DerhamOptions(degree=(2, 1, 1), bcs=(("dirichlet", "dirichlet"), None, None)), + comm=comm, + ) + Propagator.derham = derham + Propagator.domain = domain + Propagator.mass_ops = WeightedMassOperators(derham, domain) + + phi = FEECVariable(space="H1") + phi.allocate(derham=derham, domain=domain) + poisson = PoissonSolve(rho=lambda e1, e2, e3: 1.0 + 0.0 * e1) + poisson.variables.phi = phi + poisson.allocate() + assert not poisson.operator_is_singular + assert poisson._ones is None + + if __name__ == "__main__": # direction = 0 # bc_type = "dirichlet"