tests: Check get_neighbors() on all cell systems

src/python/espressomd/cellsystem.pyx

@@ -178,7 +178,21 @@ cdef class CellSystem:
 
     def get_neighbors(self, particle, distance):
         """
-        Get neighbors of a given particle up to distance
+        Get neighbors of a given particle up to a certain distance.
+
+        The choice of :ref:`cell systems <Cellsystems>` has an impact on
+        how to algorithm detects particle pairs:
+
+        * N-square: no restriction on the search distance, no double counting
+          if search distance is larger than the box size
+        * regular decomposition: the search distance is bounded by half
+          the local cell geometry
+        * hybrid decomposition: no restriction on the search distance;
+          particles in a regular decomposition cell can detect particles
+          from the same regular decomposition cell, from neighboring regular
+          decomposition cells, and from all N-square cells; particles in a
+          N-square cell can detect particles from all N-square cells and
+          from the regular decomposition cell on the same node
 
         Parameters
         ----------

testsuite/python/cellsystem_get_neighbors.py

@@ -17,16 +17,24 @@
 import unittest as ut
 import numpy as np
 import espressomd
+import espressomd.lees_edwards
 
 
-class AnalyzeDistance(ut.TestCase):
+class Test(ut.TestCase):
     system = espressomd.System(box_l=[1., 1., 1.])
     system.cell_system.skin = 0.01
     system.time_step = 0.01
     n_nodes = system.cell_system.get_state()["n_nodes"]
+    node_grid = system.cell_system.node_grid
+
+    def setUp(self):
+        self.system.cell_system.set_regular_decomposition()
+        self.system.cell_system.node_grid = self.node_grid
+        self.system.periodicity = [True, True, True]
 
     def tearDown(self):
         self.system.part.clear()
+        self.system.lees_edwards.protocol = None
 
     def get_neighbors_list(self, p1, dist):
         result = []
@@ -76,30 +84,6 @@ def test_get_neighbors_different_domains(self):
             np.testing.assert_array_equal(get_neighbors(p1, 0.1), [2])
             np.testing.assert_array_equal(get_neighbors(p2, 0.1), [1])
 
-    @ut.skipIf(n_nodes == 1, "")
-    def test_get_neighbors_greater_distance(self):
-        system = self.system
-        get_neighbors = system.cell_system.get_neighbors
-        local_box_l = np.copy(system.box_l / system.cell_system.node_grid)
-        distance = 1.001 * np.min(local_box_l)
-        pos_1 = np.array([0, 0, 0])
-        pos_2 = pos_1 + 3 * [distance]
-        p1, _ = system.part.add(pos=[pos_1, pos_2])
-        with np.testing.assert_raises_regex(Exception, r"Error in get_neighbors\(\): ERROR: pair search distance .*"):
-            get_neighbors(p1, distance)
-
-    @ut.skipIf(n_nodes in [1, 3], "Only runs for more than 1 core")
-    def test_get_neighbors_lower_distance(self):
-        system = self.system
-        # Distance needs to be lower than decomposition range
-        local_box_l = np.copy(
-            system.box_l / (2 * system.cell_system.node_grid))
-        distance = 0.9 * np.min(local_box_l)
-        pos_1 = np.array([0, 0, 0])
-        pos_2 = pos_1 + 3 * [distance]
-        system.part.add(pos=[pos_1, pos_2])
-        self.check_neighbors(distance)
-
     @ut.skipIf(n_nodes == 3, "Only runs if number of MPI cores is not 3")
     def test_get_neighbors_random_positions(self):
         system = self.system
@@ -109,6 +93,95 @@ def test_get_neighbors_random_positions(self):
             self.check_neighbors(0.2)
         system.part.clear()
 
+    def test_n_square(self):
+        """
+        Check the N-square system is supported. The search distance is not
+        bounded by the box size. No double counting occurs.
+        """
+        system = self.system
+        system.cell_system.set_n_square()
+        p1 = system.part.add(pos=[0., 0., 0.])
+        p2 = system.part.add(pos=[0., 0., 0.49])
+        self.assertEqual(system.cell_system.get_neighbors(p1, 0.5), [p2.id])
+        self.assertEqual(system.cell_system.get_neighbors(p1, 50.), [p2.id])
+
+    @ut.skipIf(n_nodes == 1, "Only runs if number of MPI cores is 2 or more")
+    def test_domain_decomposition(self):
+        """
+        Check the neighbor search distance is bounded by the regular
+        decomposition maximal range and that aperiodic systems are supported.
+        """
+        system = self.system
+        get_neighbors = system.cell_system.get_neighbors
+        local_box_l = np.copy(system.box_l / system.cell_system.node_grid)
+        max_range = np.min(local_box_l) / 2.
+        pos = np.array([0.01, 0.01, 0.01])
+        p1, p2 = system.part.add(pos=[pos, -pos])
+        self.check_neighbors(0.9999 * max_range)
+        with np.testing.assert_raises_regex(Exception, r"Error in get_neighbors\(\): ERROR: pair search distance .*"):
+            get_neighbors(p1, 1.0001 * max_range)
+
+        # change periodicity
+        pair_dist = np.linalg.norm(system.distance_vec(p1, p2))
+        self.assertEqual(get_neighbors(p1, 1.1 * pair_dist), [p2.id])
+        system.periodicity = [False, True, True]
+        self.assertEqual(get_neighbors(p1, 1.1 * pair_dist), [])
+
+    @ut.skipIf(n_nodes != 4, "Only runs if number of MPI cores is 4")
+    def test_hybrid_decomposition(self):
+        """
+        Set up colloids in a N-square cell system coupled to ions on a
+        regular decomposition cell system. The ions are set two cells
+        apart along the x-coordinate, such that they cannot see one
+        another, even when the neighbor search distance is larger than
+        the ion-ion distance.
+        """
+        system = self.system
+        system.cell_system.set_hybrid_decomposition(
+            n_square_types={1}, cutoff_regular=0.1)
+        system.cell_system.node_grid = [self.n_nodes, 1, 1]
+
+        types = {"ion": 0, "colloid": 1}
+        ions = []
+        colloids = []
+        for i in range(4):
+            pos = [0.125 + i * 0.25, 0., 0.]
+            ions.append(system.part.add(pos=pos, type=types["ion"]))
+        for i in range(4):
+            pos = [0.5, 0., 0.]
+            colloids.append(system.part.add(pos=pos, type=types["colloid"]))
+        colloid_ids = [p.id for p in colloids]
+
+        # ions can see all colloids
+        for p_ion in ions:
+            neighbors = system.cell_system.get_neighbors(p_ion, 0.5)
+            np.testing.assert_array_equal(np.sort(neighbors), colloid_ids)
+
+        # colloids can see all colloids, but only ions on the same node
+        for p_colloid in colloids:
+            neighbors = system.cell_system.get_neighbors(p_colloid, 0.5)
+            expected = list(colloid_ids)
+            expected.remove(p_colloid.id)
+            for p_ion in ions:
+                if p_ion.node == p_colloid.node:
+                    expected.insert(0, p_ion.id)
+            np.testing.assert_array_equal(np.sort(neighbors), expected)
+
+    def test_lees_edwards(self):
+        """
+        Check the Lees-Edwards position offset is taken into account
+        in the distance calculation.
+        """
+        system = self.system
+        system.lees_edwards.protocol = espressomd.lees_edwards.LinearShear(
+            initial_pos_offset=0.1, time_0=0., shear_velocity=1.0)
+        system.lees_edwards.shear_direction = 0
+        system.lees_edwards.shear_plane_normal = 2
+        p1 = system.part.add(pos=[0., 0., 0.])
+        p2 = system.part.add(pos=[0., 0., 0.99])
+        self.assertEqual(system.cell_system.get_neighbors(p1, 0.10), [])
+        self.assertEqual(system.cell_system.get_neighbors(p1, 0.15), [p2.id])
+
 
 if __name__ == "__main__":
     ut.main()