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.}])