Calculate Short Range Forces Without Ghost Force Reduction

src/core/cell_system/CellStructure.cpp

@@ -250,12 +250,14 @@ void CellStructure::set_atom_decomposition() {
 }
 
 void CellStructure::set_regular_decomposition(
-    double range, std::optional<std::pair<int, int>> fully_connected_boundary) {
+    double range, std::optional<std::pair<int, int>> fully_connected_boundary,
+    bool without_ghost_force_reduction) {
   auto &system = get_system();
   auto &local_geo = *system.local_geo;
   auto const &box_geo = *system.box_geo;
   set_particle_decomposition(std::make_unique<RegularDecomposition>(
-      ::comm_cart, range, box_geo, local_geo, fully_connected_boundary));
+      ::comm_cart, range, box_geo, local_geo, fully_connected_boundary,
+      without_ghost_force_reduction));
   m_type = CellStructureType::REGULAR;
   local_geo.set_cell_structure_type(m_type);
   system.on_cell_structure_change();

src/core/cell_system/CellStructure.hpp

@@ -558,10 +558,11 @@ struct CellStructure : public System::Leaf<CellStructure> {
    *
    * @param range Interaction range.
    * @param fully_connected_boundary neighbor cell directions for Lees-Edwards.
+   * @param without_ghost_force_reduction remove the ghost force reduction.
    */
   void set_regular_decomposition(
-      double range,
-      std::optional<std::pair<int, int>> fully_connected_boundary);
+      double range, std::optional<std::pair<int, int>> fully_connected_boundary,
+      bool without_ghost_force_reduction);
 
   /**
    * @brief Set the particle decomposition to @ref HybridDecomposition.

src/core/cell_system/HybridDecomposition.cpp

@@ -48,8 +48,9 @@ HybridDecomposition::HybridDecomposition(boost::mpi::communicator comm,
                                          LocalBox const &local_box,
                                          std::set<int> n_square_types)
     : m_comm(std::move(comm)), m_box(box_geo), m_cutoff_regular(cutoff_regular),
-      m_regular_decomposition(RegularDecomposition(
-          m_comm, cutoff_regular + skin, m_box, local_box, std::nullopt)),
+      m_regular_decomposition(
+          RegularDecomposition(m_comm, cutoff_regular + skin, m_box, local_box,
+                               std::nullopt, false)),
       m_n_square(AtomDecomposition(m_comm, m_box)),
       m_n_square_types(std::move(n_square_types)),
       m_get_global_ghost_flags(std::move(get_ghost_flags)) {

src/core/cell_system/RegularDecomposition.cpp

@@ -520,10 +520,18 @@ void RegularDecomposition::init_cell_interactions() {
 
           auto cell = &cells.at(
               get_linear_index(local_index(neighbor), ghost_cell_grid));
+
           if (ind2 > ind1) {
             red_neighbors.push_back(cell);
           } else {
-            black_neighbors.push_back(cell);
+            if (m_without_ghost_force_reduction and
+                (neighbor[2] == start[2] - 1 or neighbor[2] == end[2] or
+                 neighbor[1] == start[1] - 1 or neighbor[1] == end[1] or
+                 neighbor[0] == start[0] - 1 or neighbor[0] == end[0])) {
+              red_neighbors.push_back(cell);
+            } else {
+              black_neighbors.push_back(cell);
+            }
           }
         }
 
@@ -570,6 +578,7 @@ void assign_prefetches(GhostCommunicator &comm) {
 } // namespace
 
 GhostCommunicator RegularDecomposition::prepare_comm() {
+
   int dir, lr, i, cnt, n_comm_cells[3];
   Utils::Vector3i lc{}, hc{}, done{};
 
@@ -670,10 +679,11 @@ GhostCommunicator RegularDecomposition::prepare_comm() {
 RegularDecomposition::RegularDecomposition(
     boost::mpi::communicator comm, double range, BoxGeometry const &box_geo,
     LocalBox const &local_geo,
-    std::optional<std::pair<int, int>> fully_connected)
+    std::optional<std::pair<int, int>> fully_connected,
+    bool without_ghost_force_reduction)
     : m_comm(std::move(comm)), m_box(box_geo), m_local_box(local_geo),
-      m_fully_connected_boundary(std::move(fully_connected)) {
-
+      m_fully_connected_boundary(std::move(fully_connected)),
+      m_without_ghost_force_reduction(without_ghost_force_reduction) {
   /* set up new regular decomposition cell structure */
   create_cell_grid(range);
 
@@ -685,11 +695,13 @@ RegularDecomposition::RegularDecomposition(
 
   /* create communicators */
   m_exchange_ghosts_comm = prepare_comm();
-  m_collect_ghost_force_comm = prepare_comm();
-
-  /* collect forces has to be done in reverted order! */
-  revert_comm_order(m_collect_ghost_force_comm);
-
   assign_prefetches(m_exchange_ghosts_comm);
-  assign_prefetches(m_collect_ghost_force_comm);
+  if (m_without_ghost_force_reduction) {
+    m_collect_ghost_force_comm = GhostCommunicator{};
+  } else {
+    m_collect_ghost_force_comm = prepare_comm();
+    /* collect forces has to be done in reverted order! */
+    revert_comm_order(m_collect_ghost_force_comm);
+    assign_prefetches(m_collect_ghost_force_comm);
+  }
 }

src/core/cell_system/RegularDecomposition.hpp

@@ -85,12 +85,18 @@ struct RegularDecomposition : public ParticleDecomposition {
   std::vector<Cell *> m_ghost_cells;
   GhostCommunicator m_exchange_ghosts_comm;
   GhostCommunicator m_collect_ghost_force_comm;
+  bool m_without_ghost_force_reduction;
 
 public:
-  RegularDecomposition(boost::mpi::communicator comm, double range,
-                       BoxGeometry const &box_geo, LocalBox const &local_geo,
-                       std::optional<std::pair<int, int>> fully_connected);
-
+  RegularDecomposition(
+      boost::mpi::communicator comm, double range, BoxGeometry const &box_geo,
+      LocalBox const &local_geo,
+      std::optional<std::pair<int, int>> fully_connected_boundary,
+      bool without_ghost_force_reduction = false);
+
+  bool get_without_ghost_force_reduction() const noexcept {
+    return m_without_ghost_force_reduction;
+  }
   GhostCommunicator const &exchange_ghosts_comm() const override {
     return m_exchange_ghosts_comm;
   }

src/core/cells.cpp

@@ -201,4 +201,4 @@ std::vector<NeighborPIDs> get_neighbor_pids(System::System const &system) {
     kernel(p, get_interacting_neighbors(system, p));
   }
   return ret;
-}
+}

src/core/system/System.cpp

@@ -162,9 +162,11 @@ void System::set_cell_structure_topology(CellStructureType topology) {
               std::as_const(*cell_structure).decomposition());
       cell_structure->set_regular_decomposition(
           get_interaction_range(),
-          old_regular_decomposition.fully_connected_boundary());
+          old_regular_decomposition.fully_connected_boundary(),
+          old_regular_decomposition.get_without_ghost_force_reduction());
     } else { // prev. decomposition is not a regular decomposition
-      cell_structure->set_regular_decomposition(get_interaction_range(), {});
+      cell_structure->set_regular_decomposition(get_interaction_range(), {},
+                                                false);
     }
   } else if (topology == CellStructureType::NSQUARE) {
     cell_structure->set_atom_decomposition();

src/script_interface/cell_system/CellSystem.cpp

@@ -162,6 +162,15 @@ CellSystem::CellSystem() {
        [this]() { return get_system().bonded_ias->maximal_cutoff(); }},
       {"interaction_range", AutoParameter::read_only,
        [this]() { return get_system().get_interaction_range(); }},
+      {"without_ghost_force_reduction", AutoParameter::read_only,
+       [this]() {
+         if (get_cell_structure().decomposition_type() !=
+             CellStructureType::REGULAR) {
+           return Variant{none};
+         }
+         auto const rd = get_regular_decomposition();
+         return Variant{rd.get_without_ghost_force_reduction()};
+       }},
   });
 }
 
@@ -292,6 +301,14 @@ void CellSystem::initialize(CellStructureType const &cs_type,
   auto &system = get_system();
   m_cell_structure->use_verlet_list = verlet;
   if (cs_type == CellStructureType::HYBRID) {
+    auto without_ghost_force_reduction =
+        get_value_or<bool>(params, "without_ghost_force_reduction", false);
+    context()->parallel_try_catch([without_ghost_force_reduction]() {
+      if (without_ghost_force_reduction) {
+        throw std::invalid_argument("Parameter 'without_ghost_force_reduction' "
+                                    "is not allowed for hybrid decomposition");
+      }
+    });
     auto const cutoff_regular = get_value<double>(params, "cutoff_regular");
     auto const ns_types =
         get_value_or<std::vector<int>>(params, "n_square_types", {});
@@ -308,10 +325,14 @@ void CellSystem::initialize(CellStructureType const &cs_type,
                      coord(boost::get<std::string>(variant.at("direction")))}};
       });
     }
-    context()->parallel_try_catch([this, &fcb_pair]() {
-      m_cell_structure->set_regular_decomposition(
-          get_system().get_interaction_range(), fcb_pair);
-    });
+    auto const without_ghost_force_reduction =
+        get_value_or<bool>(params, "without_ghost_force_reduction", false);
+    context()->parallel_try_catch(
+        [this, &fcb_pair, without_ghost_force_reduction]() {
+          m_cell_structure->set_regular_decomposition(
+              get_system().get_interaction_range(), fcb_pair,
+              without_ghost_force_reduction);
+        });
   } else {
     system.set_cell_structure_topology(cs_type);
   }

testsuite/python/cell_system.py

@@ -83,14 +83,17 @@ def test_exceptions(self):
             system.cell_system.node_grid = [1, 2, self.n_nodes]
         np.testing.assert_array_equal(np.copy(system.cell_system.node_grid),
                                       np.copy(node_grid))
+        with self.assertRaisesRegex(ValueError, "Parameter 'without_ghost_force_reduction' is not allowed for hybrid decomposition"):
+            system.cell_system.set_hybrid_decomposition(
+                n_square_types={1}, cutoff_regular=0, without_ghost_force_reduction=True)
 
     def test_node_grid_regular(self):
         self.system.cell_system.set_regular_decomposition()
         self.check_node_grid()
 
     def test_node_grid_hybrid(self):
         self.system.cell_system.set_hybrid_decomposition(
-            n_square_types={1}, cutoff_regular=0)
+            n_square_types={1}, cutoff_regular=0, without_ghost_force_reduction=False)
         self.check_node_grid()
 
 

testsuite/python/lj.py

@@ -72,6 +72,18 @@ def test_dd_vl(self):
         self.system.integrator.run(recalc_forces=True, steps=0)
         self.check()
 
+    def test_dd_with_ghost_force_reduction(self):
+        self.system.cell_system.set_regular_decomposition(
+            without_ghost_force_reduction=False)
+        self.system.integrator.run(recalc_forces=True, steps=0)
+        self.check()
+
+    def test_dd_without_ghost_force_reduction(self):
+        self.system.cell_system.set_regular_decomposition(
+            without_ghost_force_reduction=True)
+        self.system.integrator.run(recalc_forces=True, steps=0)
+        self.check()
+
 
 if __name__ == '__main__':
     ut.main()