From 1d2e28e422e4665ba753c8f869d833197d2acad3 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Tue, 13 Oct 2026 08:40:12 +0200 Subject: [PATCH] draw_markers: send the markers to their processes also without sorting boxes (#764) On more than one process, draw_markers MPI-sorted the markers only if sorting boxes were set up. Without boxes, the markers drawn on the whole domain stayed on the wrong process until the first push, and accumulations before it (e.g. the initial Poisson solve of VlasovAmpereOneSpecies) wrote outside the local stencil arrays (segfault on 2 ranks). Now draw_markers always MPI-sorts on several processes; sorting into boxes is unchanged. On one process the behaviour is unchanged. New MPI test test_draw_without_boxes_on_rank_domain (pseudo_random, sobol_standard, sobol_antithetic; sort=True/False): fails on 2 ranks without the fix. Closes #764. Co-Authored-By: Claude Opus 5.5 --- src/struphy/pic/base.py | 19 +++++--- src/struphy/pic/tests/test_draw_parallel.py | 50 +++++++++++++++++++++ 2 files changed, 63 insertions(+), 6 deletions(-) diff --git a/src/struphy/pic/base.py b/src/struphy/pic/base.py index 5e1d217b2..c5c629389 100644 --- a/src/struphy/pic/base.py +++ b/src/struphy/pic/base.py @@ -1146,12 +1146,14 @@ def draw_markers( All needed parameters can be set in the parameter file, see :ref:`params_yml`. - An initial sorting will be performed if sort is given as True (default) and sorting_params were given to the init. + On several processes, the markers are always sent to the process of their domain after drawing + (:meth:`mpi_sort_markers`). They are also sorted into boxes if sort is True (default) and sorting + boxes were set up. Parameters ---------- sort : Bool - Wether to sort the particules in boxes after initial drawing (only if sorting params were passed) + Whether to sort the markers into boxes after initial drawing (only if sorting boxes were set up). """ import scipy.special as sp @@ -1405,10 +1407,15 @@ def draw_markers( assert xp.all(~self.holes[:n_mks_load_loc]) assert xp.all(self.holes[n_mks_load_loc:]) - if self._initialized_sorting and sort: - logger.info("\nSorting the markers after initial draw") - if self.mpi_comm is not None: - self.mpi_sort_markers() + # send the markers to the process of their domain: needed on several processes whether or not there are + # sorting boxes, since accumulations before the first push (e.g. an initial Poisson solve) assume it (#764) + sort_boxes = self._initialized_sorting and sort + if self.mpi_comm is not None and (self.mpi_size > 1 or sort_boxes): + logger.info("\nSending the markers to their processes after initial draw") + self.mpi_sort_markers() + + if sort_boxes: + logger.info("\nSorting the markers into boxes after initial draw") self.do_sort() logger.info("Done.") diff --git a/src/struphy/pic/tests/test_draw_parallel.py b/src/struphy/pic/tests/test_draw_parallel.py index ff18a3321..09020c588 100644 --- a/src/struphy/pic/tests/test_draw_parallel.py +++ b/src/struphy/pic/tests/test_draw_parallel.py @@ -177,6 +177,56 @@ def test_marker_ids(loading, Np=1000): assert n_at_origin == 0 +@pytest.mark.parametrize("loading", ["pseudo_random", "sobol_standard", "sobol_antithetic"]) +@pytest.mark.parametrize("sort", [True, False]) +def test_draw_without_boxes_on_rank_domain(loading, sort, Np=1000): + """Without sorting boxes, every marker is on the domain of its process right after drawing (#764), + for ``sort=True`` and ``sort=False`` (which only concerns the boxes), and no marker is lost.""" + + import cunumpy as xp + from maybempi import MPI + + from struphy import LoadingParameters, domains + from struphy.feec.psydac_derham import Derham + from struphy.io.options import DerhamOptions + from struphy.pic.particles import Particles6D + from struphy.topology.grids import TensorProductGrid + + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=[8, 5, 3]) + derham = Derham(grid, DerhamOptions(degree=[1, 1, 1]), comm=comm) + domain_decomp = (derham.domain_array, derham.domain_decomposition.nprocs) + + loading_params = LoadingParameters( + Np=Np, + seed=1234, + loading=loading, + moments=(0.0, 0.0, 0.0, 1.0, 1.0, 1.0), + ) + + # no sorting_params: no sorting boxes + particles = Particles6D( + comm_world=comm, + domain_decomp=domain_decomp, + loading_params=loading_params, + domain=domains.Cuboid(), + ) + assert particles.sorting_boxes is None + + particles.draw_markers(sort=sort) + + positions = particles.markers_wo_holes[:, :3] + on_rank = xp.all(particles._is_on_domain_of_rank(positions, particles.mpi_rank), axis=1) + n_off_rank = int(xp.sum(~on_rank)) + n_markers = positions.shape[0] + if particles.mpi_size > 1: + n_markers = particles.mpi_comm.allreduce(n_markers) + + assert n_off_rank == 0, f"rank {particles.mpi_rank}: {n_off_rank} markers off the process domain after drawing" + assert n_markers == Np + + if __name__ == "__main__": # test_draw([8, 9, 10], [2, 3, 4], [False, False, True], ['Cuboid', { # 'l1': 1., 'r1': 2., 'l2': 10., 'r2': 20., 'l3': 100., 'r3': 200.}])