From 4d4b51f35db87fee555282df8364bb38ed165db7 Mon Sep 17 00:00:00 2001 From: Max Date: Sat, 10 Oct 2026 08:18:35 +0200 Subject: [PATCH 1/2] Diagonal scaling: include band entries that wrap onto the diagonal In a periodic direction with fewer points than the band width (2p + 1 > n, e.g. a single element in a 1D run), StencilMatrix.diagonal() returns only the entry at offset 0, while the entries at offsets k = m n also lie on the main diagonal. The diagonal D of the mass matrix was therefore wrong, and MassMatrixPreconditioner with diagonal_scaling=True (the default since #739) was no longer exact on the logical cube; in VlasovAmpereCoupling this broke energy conservation in 1D runs. Co-Authored-By: Claude Opus 5.5 --- src/struphy/feec/preconditioner.py | 30 +++++++++++++++-- src/struphy/feec/tests/test_mass_matrices.py | 34 ++++++++++++++++++++ 2 files changed, 62 insertions(+), 2 deletions(-) diff --git a/src/struphy/feec/preconditioner.py b/src/struphy/feec/preconditioner.py index d64838dc2..a3a67c70b 100644 --- a/src/struphy/feec/preconditioner.py +++ b/src/struphy/feec/preconditioner.py @@ -318,7 +318,8 @@ def weight_1d(c: int, d: int) -> Callable | xp.ndarray: return _kronecker_approximation(mass_operator, bc, weight_1d) def _core_diagonal(self): - return self._mass_operator.matrix.diagonal() + # the matrix itself: _local_diagonal also adds the band entries that wrap onto the diagonal + return self._mass_operator.matrix @property def dim_reduce(self) -> int | None: @@ -765,7 +766,7 @@ def _local_diagonal(op: LinearOperator) -> list[xp.ndarray]: if isinstance(op, StencilDiagonalMatrix): return [op._data] if isinstance(op, StencilMatrix): - return [op.diagonal()._data] + return [_stencil_local_diagonal(op)] if isinstance(op, KroneckerStencilMatrix): # outer product of the diagonals of the (process-local) factors on the local rows diags = [] @@ -788,6 +789,31 @@ def _local_diagonal(op: LinearOperator) -> list[xp.ndarray]: raise NotImplementedError(f"Diagonal of {type(op).__name__} is not supported.") +def _stencil_local_diagonal(M: StencilMatrix) -> xp.ndarray: + """ + Local diagonal of a StencilMatrix, including the band entries that wrap onto the main diagonal. + + In a periodic direction with fewer points than the band width (``2p + 1 > n``, e.g. a single + element), the band entries at the offsets ``k = m n`` (``m != 0``) are on the main diagonal too; + :meth:`StencilMatrix.diagonal` returns only the entry at ``k = 0``. + """ + V = M.codomain + assert all(m == 1 for m in V.shifts) and all(m == 1 for m in M.domain.shifts), "Only shifts 1 are supported." + offsets = [ + [k for k in range(-p, p + 1) if k % n == 0] if periodic else [0] + for n, p, periodic in zip(V.npts, M.pads, V.periods) + ] + diag = M.diagonal()._data + if all(o == [0] for o in offsets): + return diag + idx = M._get_diagonal_indices() + rows, cols = idx[: V.ndim], idx[V.ndim :] + for ks in itertools.product(*offsets): + if any(ks): + diag = diag + M._data[rows + tuple(j + k for j, k in zip(cols, ks))] + return diag + + def _diagonal_operator(space: VectorSpace, diags: list[xp.ndarray]) -> LinearOperator: """Diagonal operator on ``space`` (Stencil- or BlockVectorSpace) with the given local diagonals.""" if isinstance(space, BlockVectorSpace): diff --git a/src/struphy/feec/tests/test_mass_matrices.py b/src/struphy/feec/tests/test_mass_matrices.py index 23af76b82..799646edf 100644 --- a/src/struphy/feec/tests/test_mass_matrices.py +++ b/src/struphy/feec/tests/test_mass_matrices.py @@ -1742,6 +1742,40 @@ def test_transpose_and_copy(): assert xp.allclose(A.T.toarray(), -A_arr, atol=1e-14) +@pytest.mark.parametrize( + "num_elements, degree", [((16, 1, 1), (3, 1, 1)), ((16, 2, 1), (2, 2, 1)), ((8, 6, 4), (2, 2, 2))] +) +@pytest.mark.parametrize("space", ["M0", "M1", "M2", "M3"]) +def test_mass_preconditioner_diagonal_scaling_periodic_exact(num_elements, degree, space): + """On the periodic logical cube the mass matrices are Kronecker products, so MassMatrixPreconditioner with + diagonal scaling inverts them exactly; also with fewer points than the band width in a direction (one element), + where band entries wrap onto the main diagonal.""" + + import numpy as np + from maybempi import MPI + + from struphy import domains + from struphy.feec.mass import WeightedMassOperators + from struphy.feec.preconditioner import MassMatrixPreconditioner + from struphy.feec.psydac_derham import Derham + from struphy.io.options import DerhamOptions + from struphy.topology.grids import TensorProductGrid + + domain = domains.Cuboid(r1=12.56) + derham = Derham(TensorProductGrid(num_elements=num_elements), DerhamOptions(degree=degree), comm=MPI.COMM_WORLD) + M = getattr(WeightedMassOperators(derham, domain), space) + + rng = np.random.default_rng(0) + v = M.domain.zeros() + for block in v.blocks if hasattr(v, "blocks") else [v]: + block._data[:] = rng.random(block._data.shape) + v.update_ghost_regions() + + w = MassMatrixPreconditioner(M, diagonal_scaling=True).dot(M.dot(v)) + err = ((w - v).inner(w - v) / v.inner(v)) ** 0.5 + assert err < 1e-12, f"{space}, {num_elements=}, {degree=}: |P M v - v| / |v| = {err}" + + if __name__ == "__main__": # test_mass( # num_elements=(32, 32, 32), From a0f1cf35558b5c985728fbad4b60915e85aa9bed Mon Sep 17 00:00:00 2001 From: Max Date: Sat, 10 Oct 2026 08:23:55 +0200 Subject: [PATCH 2/2] VlasovAmpereCoupling: update the ghost regions of the mid-point field before the push The pusher kernel reads the ghost regions of e_sum directly. They were only correct if the Schur solver happened to return its solution with valid ghost regions, which is no longer the case with the diagonal scaling of MassMatrixPreconditioner (default since #739): its last step is a StencilDiagonalMatrix.dot that writes only the owned entries. In 1D runs (one element in e2, e3) every evaluation reads ghost values, and the total energy was no longer conserved (drift ~1e-4 instead of ~1e-12). The other couplings already update the ghost regions before pushing. Co-Authored-By: Claude Opus 5.5 --- src/struphy/propagators/vlasov_ampere_coupling.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/struphy/propagators/vlasov_ampere_coupling.py b/src/struphy/propagators/vlasov_ampere_coupling.py index 0c8c76414..77a129ae0 100644 --- a/src/struphy/propagators/vlasov_ampere_coupling.py +++ b/src/struphy/propagators/vlasov_ampere_coupling.py @@ -214,6 +214,8 @@ def __call__(self, dt): self._e_sum += self.variables.e.spline.vector self._e_sum += self._e_tmp self._e_sum *= 0.5 + # the pusher kernel reads the ghost regions directly (the solver output may leave them stale) + self._e_sum.update_ghost_regions() # Update velocities self._pusher(dt)