Draft: Add flag for 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,6 +520,7 @@ 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 {
@@ -670,10 +671,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);
 

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()};
+       }},
   });
 }
 
@@ -299,6 +308,10 @@ void CellSystem::initialize(CellStructureType const &cs_type,
     m_cell_structure->set_hybrid_decomposition(cutoff_regular, n_square_types);
   } else if (cs_type == CellStructureType::REGULAR) {
     std::optional<std::pair<int, int>> fcb_pair = std::nullopt;
+    if (get_value_or(params, "without_ghost_force_reduction", false)) {
+      throw std::invalid_argument("Parameter 'without_ghost_force_reduction' "
+                                  "is not allowed for hybrid decomposition");
+    }
     if (params.contains("fully_connected_boundary") and
         not is_none(params.at("fully_connected_boundary"))) {
       auto const variant =
@@ -308,10 +321,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);
   }