Make VS relative thread-safe (#5212)

Author: RudolfWeeber

src/core/cell_system/CellStructure.hpp

@@ -364,9 +364,10 @@ class CellStructure : public System::Leaf<CellStructure> {
    * @brief Run a kernel on all local particles.
    * The kernel is assumed to be thread-safe.
    */
-  void for_each_local_particle(ParticleUnaryOp &&f) const {
+  void for_each_local_particle(ParticleUnaryOp &&f,
+                               bool parallel = true) const {
 #ifdef ESPRESSO_SHARED_MEMORY_PARALLELISM
-    if (use_parallel_for_each_local_particle()) {
+    if (parallel and use_parallel_for_each_local_particle()) {
       parallel_for_each_particle_impl(decomposition().local_cells(), f);
       return;
     }

src/core/virtual_sites/relative.cpp

@@ -168,23 +168,26 @@ void vs_relative_back_transfer_forces_and_torques(
   init_forces_ghosts(cell_structure);
 
   // Iterate over all the particles in the local cells
-  cell_structure.for_each_local_particle([&](Particle &p) {
-    if (!is_vs(p))
-      return;
-
-    auto *p_ref_ptr = get_reference_particle(cell_structure, p);
-    assert(p_ref_ptr != nullptr);
-
-    auto &p_ref = *p_ref_ptr;
-    if (is_vs_relative_trans(p)) {
-      p_ref.force() += p.force();
-      p_ref.torque() += vector_product(connection_vector(p_ref, p), p.force());
-    }
-
-    if (is_vs_rot(p)) {
-      p_ref.torque() += p.torque();
-    }
-  });
+  cell_structure.for_each_local_particle(
+      [&](Particle &p) {
+        if (!is_vs(p))
+          return;
+
+        auto *p_ref_ptr = get_reference_particle(cell_structure, p);
+        assert(p_ref_ptr != nullptr);
+
+        auto &p_ref = *p_ref_ptr;
+        if (is_vs_relative_trans(p)) {
+          p_ref.force() += p.force();
+          p_ref.torque() +=
+              vector_product(connection_vector(p_ref, p), p.force());
+        }
+
+        if (is_vs_rot(p)) {
+          p_ref.torque() += p.torque();
+        }
+      },
+      /* parallel */ false);
 }
 
 // Rigid body contribution to scalar pressure and pressure tensor

testsuite/python/virtual_sites_relative.py

@@ -202,14 +202,19 @@ def test_pos_vel_forces(self):
         system.min_global_cut = 0.23
         self.assertEqual(system.min_global_cut, 0.23)
 
-        # Place central particle + 3 vs
+        # Place central particle + N virtual sites
         p1 = system.part.add(rotation=3 * [True], pos=(0.5, 0.5, 0.5), id=1,
                              quat=(1, 0, 0, 0), omega_lab=(1, 2, 3))
-        pos2 = (0.5, 0.4, 0.5)
-        pos3 = (0.3, 0.5, 0.4)
-        pos4 = (0.5, 0.5, 0.5)
-        for pos in (pos2, pos3, pos4):
-            p = system.part.add(rotation=3 * [True], pos=pos)
+
+        # Number of virtual sites to create
+        N = 100 
+        # Generate N random positions within 1.2 of central particle in each coordinate
+        np.random.seed(42)
+        vs_positions = p1.pos + np.random.uniform(-0.15, 0.15, (N, 3))
+
+        # Create virtual sites at random positions
+        sites = system.part.add(rotation=[3 * [True]] * N, pos=vs_positions)
+        for p in sites: 
             p.vs_auto_relate_to(p1)
             # Was the particle made virtual
             self.assertTrue(p.is_virtual())
@@ -226,39 +231,39 @@ def test_pos_vel_forces(self):
         p1.v = (0.45, 0.14, 0.447)
         p1.omega_lab = (0.45, 0.14, 0.447)
         system.integrator.run(0, recalc_forces=True)
-        for p in system.part:
-            if p.id != p1.id:
-                self.verify_vs(p)
+        for p in sites:
+            self.verify_vs(p)
 
         # Check if still true, when non-virtual particle has rotated and a
         # linear motion
         p1.omega_lab = [-5., 3., 8.4]
         system.integrator.run(10)
-        for p in system.part:
-            if p.id != p1.id:
-                self.verify_vs(p)
+        for p in sites:
+            self.verify_vs(p)
 
         if espressomd.has_features("EXTERNAL_FORCES"):
             # Test transfer of forces accumulating on virtual sites
             # to central particle
-            f2 = np.array((3, 4, 5))
-            f3 = np.array((-4, 5, 6))
-            # Add forces to vs
-            p2, p3 = system.part.by_ids([2, 3])
-            p2.ext_force = f2
-            p3.ext_force = f3
+            # Generate random forces for all N virtual sites
+            sites.ext_force = np.random.uniform(-5, 5, (N, 3))
+
             system.integrator.run(0)
             # get force/torques on non-vs
             f = p1.f
             t = p1.torque_lab
 
-            # Expected force = sum of the forces on the vs
-            self.assertAlmostEqual(np.linalg.norm(f - f2 - f3), 0., delta=1E-6)
+            # Expected force = sum of all forces on the vs
+            f_exp = np.sum(sites.ext_force, axis=0)
+            print()
+            print(f"{f=} {p1.f=}, {f_exp=}")
+            self.assertAlmostEqual(np.linalg.norm(f - f_exp), 0., delta=1E-6)
 
             # Expected torque
             # Radial components of forces on a rigid body add to the torque
-            t_exp = np.cross(system.distance_vec(p1, p2), f2)
-            t_exp += np.cross(system.distance_vec(p1, p3), f3)
+            t_exp = np.zeros(3)
+            for vs_p in sites:
+                t_exp += np.cross(system.distance_vec(p1,
+                                  vs_p), vs_p.ext_force)
             # Check
             self.assertAlmostEqual(np.linalg.norm(t_exp - t), 0., delta=1E-6)