galilei: fix out-of-bounds access for duplicate ComFixed types (bug-sweep #27)

src/core/galilei/ComFixed.hpp

@@ -72,9 +72,15 @@ class ComFixed {
   void set_fixed_types(std::vector<int> const &p_types) {
     m_type_index.clear();
 
-    int i = 0;
+    /* Assign dense indices 0..n_distinct-1. Using try_emplace with the
+     * current map size as the candidate index makes duplicate types collapse
+     * onto a single index instead of inflating the counter, so every stored
+     * index stays < m_type_index.size(). Consumers (local_type_forces,
+     * local_type_masses, apply) size their vectors by m_type_index.size() and
+     * index by it->second, so this invariant is required to avoid an
+     * out-of-bounds vector access. */
     for (auto const &p_type : p_types) {
-      m_type_index[p_type] = i++;
+      m_type_index.try_emplace(p_type, static_cast<int>(m_type_index.size()));
     }
   }
 

testsuite/python/comfixed.py

@@ -28,45 +28,67 @@ class ComFixed(ut.TestCase):
 
     np.random.seed(seed=42)
 
-    def com(self, parts):
-        return np.average(parts.pos, axis=0, weights=parts.mass)
+    system = espressomd.System(box_l=[10., 10., 10.])
+    system.time_step = 0.01
+    system.cell_system.skin = 0.4
 
-    def test(self):
-        dt = 0.01
-        skin = 0.4
+    def setUp(self):
+        self.system.thermostat.set_langevin(kT=1., gamma=0.01, seed=41)
 
-        system = espressomd.System(box_l=[10., 10., 10.])
-        system.time_step = dt
-        system.cell_system.skin = skin
+    def tearDown(self):
+        self.system.comfixed.types = []
+        self.system.thermostat.turn_off()
+        self.system.part.clear()
 
-        system.thermostat.set_langevin(kT=1., gamma=0.01, seed=41)
+    def com(self, parts):
+        return np.average(parts.pos, axis=0, weights=parts.mass)
+
+    def check_com_conserved(self, fixed_types, particle_types):
+        """Set up a system whose particles carry the given ``particle_types``,
+        fix the centre of mass of ``fixed_types`` and assert that the COM of
+        the fixed particles does not drift during integration."""
+        system = self.system
 
         for i in range(100):
             r = [0.5, 1., 1.] * system.box_l * np.random.random(3)
             v = 3 * [0.]
             # Make sure that id and type gaps work correctly
-            system.part.add(id=2 * i, pos=r, v=v, type=2 * (i % 2))
-        partcls = system.part.all()
+            system.part.add(id=2 * i, pos=r, v=v,
+                            type=particle_types[i % len(particle_types)])
 
         if espressomd.has_features(["MASS"]):
             # Avoid masses too small for the time step
-            partcls.mass = 2. * (0.1 + np.random.random(100))
-
-        com_0 = self.com(partcls)
+            system.part.all().mass = 2. * (0.1 + np.random.random(100))
 
-        system.comfixed.types = [0, 2]
+        # COM of the subset of particles whose type is fixed.
+        distinct_fixed = set(fixed_types)
+        fixed_parts = system.part.select(lambda p: p.type in distinct_fixed)
+        com_0 = self.com(fixed_parts)
 
-        # Interface check
-        self.assertEqual(system.comfixed.types, [2, 0])
-
-        for i in range(10):
-            com_i = self.com(partcls)
+        system.comfixed.types = fixed_types
 
+        for _ in range(10):
+            com_i = self.com(fixed_parts)
             for j in range(3):
                 self.assertAlmostEqual(com_0[j], com_i[j], places=10)
-
             system.integrator.run(100)
 
+    def test(self):
+        self.check_com_conserved(fixed_types=[0, 2], particle_types=[0, 2])
+
+        # Interface check
+        self.system.comfixed.types = [0, 2]
+        self.assertEqual(self.system.comfixed.types, [2, 0])
+
+    def test_duplicate_types(self):
+        """Regression test for bug #27: a duplicated entry in the ComFixed
+        types list used to leave the internal type->index map with a stored
+        index >= map size, producing an out-of-bounds vector access in
+        ComFixed::apply() during integration. The duplicate must be handled
+        safely (deduplicated): fixing the COM of types ``[1, 1]`` must behave
+        exactly like fixing the COM of type ``1`` and conserve it."""
+        self.check_com_conserved(fixed_types=[1, 1], particle_types=[1, 3])
+
 
 if __name__ == "__main__":
     ut.main()