diff --git a/runtime/domain/block_forest.cpp b/runtime/domain/block_forest.cpp
index 81a416641b9c9329a5c7ebadd4cf4565ef923e30..1ce021efd1f3814915e990b77ac0ef1f8e4f7486 100644
--- a/runtime/domain/block_forest.cpp
+++ b/runtime/domain/block_forest.cpp
@@ -1,3 +1,4 @@
+#include <map>
 #include <mpi.h>
 #include <vector>
 //---
@@ -21,6 +22,8 @@
 namespace pairs {
 
 void BlockForest::updateNeighborhood() {
+    std::map<int, std::vector<walberla::math::AABB>> neighborhood;
+    std::map<int, std::vector<walberla::BlockID>> blocks_pushed;
     auto me = walberla::mpi::MPIManager::instance()->rank();
     this->nranks = 0;
     this->total_aabbs = 0;
diff --git a/runtime/domain/block_forest.hpp b/runtime/domain/block_forest.hpp
index 673369d6bb925317de43717c5bc73cd5e8c3fa9c..acfa0c5a463489d4f042dbb5dde5df88a9c42f18 100644
--- a/runtime/domain/block_forest.hpp
+++ b/runtime/domain/block_forest.hpp
@@ -1,5 +1,3 @@
-#include <map>
-//---
 #include <blockforest/BlockForest.h>
 #include <blockforest/Initialization.h>
 #include <blockforest/loadbalancing/DynamicCurve.h>
@@ -27,8 +25,6 @@ private:
     PairsSimulation *ps;
     std::shared_ptr<walberla::BlockForest> forest;
     std::shared_ptr<walberla::blockforest::InfoCollection> info;
-    std::map<int, std::vector<walberla::math::AABB>> neighborhood;
-    std::map<int, std::vector<walberla::BlockID>> blocks_pushed;
     std::vector<int> ranks;
     std::vector<int> naabbs;
     std::vector<int> aabb_offsets;
diff --git a/src/pairs/sim/comm.py b/src/pairs/sim/comm.py
index 08fb1438ae62ce34d57b2f76435f33db304e6780..fbd8de778727b44d57aceab40ca5e296036bf01c 100644
--- a/src/pairs/sim/comm.py
+++ b/src/pairs/sim/comm.py
@@ -93,7 +93,7 @@ class Comm:
             CommunicateData(self, step, prop_list)
             UnpackGhostParticles(self, step, prop_list)
 
-            step_nrecv = sum([self.nrecv[j] for j in self.dom_part.step_indexes(step)])
+            step_nrecv = self.dom_part.reduce_sum_step_indexes(step, self.nrecv)
             Assign(self.sim, self.sim.nghost, self.sim.nghost + step_nrecv)
 
     @pairs_inline
@@ -306,9 +306,9 @@ class PackGhostParticles(Lowerable):
         send_mult = self.comm.send_mult
         self.sim.module_name(f"pack_ghost_particles{self.step}_" + "_".join([str(p.id()) for p in self.prop_list]))
 
-        step_indexes = self.comm.dom_part.step_indexes(self.step)
-        start = self.comm.send_offsets[step_indexes[0]]
-        for i in For(self.sim, start, ScalarOp.inline(start + sum([self.comm.nsend[j] for j in step_indexes]))):
+        start = self.comm.send_offsets[self.comm.dom_part.first_step_index(self.step)]
+        end = ScalarOp.inline(start + self.comm.dom_part.reduce_sum_step_indexes(self.step, self.comm.nsend))
+        for i in For(self.sim, start, end):
             p_offset = 0
             m = send_map[i]
             for p in self.prop_list:
@@ -347,9 +347,9 @@ class UnpackGhostParticles(Lowerable):
         recv_buffer.set_stride(1, self.get_elems_per_particle())
         self.sim.module_name(f"unpack_ghost_particles{self.step}_" + "_".join([str(p.id()) for p in self.prop_list]))
 
-        step_indexes = self.comm.dom_part.step_indexes(self.step)
-        start = self.comm.recv_offsets[step_indexes[0]]
-        for i in For(self.sim, start, ScalarOp.inline(start + sum([self.comm.nrecv[j] for j in step_indexes]))):
+        start = self.comm.recv_offsets[self.comm.dom_part.first_step_index(self.step)]
+        end = ScalarOp.inline(start + self.comm.dom_part.reduce_sum_step_indexes(self.step, self.comm.nrecv))
+        for i in For(self.sim, start, end):
             p_offset = 0
             for p in self.prop_list:
                 if not Types.is_scalar(p.type()):
@@ -422,9 +422,7 @@ class UnpackAllGhostParticles(Lowerable):
         recv_buffer.set_stride(1, self.get_elems_per_particle())
         self.sim.module_name(f"unpack_all_ghost_particles" + "_".join([str(p.id()) for p in self.prop_list]))
 
-        nrecv_size = sum([len(dom_part.step_indexes(s)) for s in range(dom_part.number_of_steps())])
-        nrecv_all = sum([self.comm.nrecv[j] for j in range(nrecv_size)])
-
+        nrecv_all = self.comm.dom_part.reduce_sum_all_steps(self.comm.nrecv)
         for i in For(self.sim, 0, nrecv_all):
             p_offset = 0
             for p in self.prop_list:
@@ -516,7 +514,8 @@ class ChangeSizeAfterExchange(Lowerable):
     def lower(self):
         self.sim.module_name(f"change_size_after_exchange{self.step}")
         self.sim.check_resize(self.sim.particle_capacity, self.sim.nlocal)
-        Assign(self.sim, self.sim.nlocal, self.sim.nlocal + sum([self.comm.nrecv[j] for j in self.comm.dom_part.step_indexes(self.step)]))
+        nrecv = self.comm.dom_part.reduce_sum_step_indexes(self.step, self.comm.nrecv)
+        Assign(self.sim, self.sim.nlocal, self.sim.nlocal + nrecv)
 
 
 class PackContactHistoryData(Lowerable):
@@ -611,7 +610,6 @@ class UnpackContactHistoryData(Lowerable):
         contact_used = self.sim._contact_history.contact_used
         self.sim.module_name(f"unpack_contact_history{self.step}")
 
-        step_indexes = self.comm.dom_part.step_indexes(self.step)
         nelems_per_contact = sum([Types.number_of_elements(self.sim, cp.type()) \
                                   for cp in self.sim.contact_properties]) + 1
 
diff --git a/src/pairs/sim/domain_partitioning.py b/src/pairs/sim/domain_partitioning.py
index 83875af3847a2a61d787f910a270435835f6051a..6cf3dbf0bf8f0ed7d5518ca6f92f2045d3bd6a16 100644
--- a/src/pairs/sim/domain_partitioning.py
+++ b/src/pairs/sim/domain_partitioning.py
@@ -29,6 +29,16 @@ class DimensionRanges:
     def step_indexes(self, step):
         return [step * 2 + 0, step * 2 + 1]
 
+    def first_step_index(self, step):
+        return self.step_indexes(step)[0]
+
+    def reduce_sum_all_steps(self, array):
+        total_size = sum([len(self.step_indexes(s)) for s in range(self.number_of_steps())])
+        return sum([array[i] for i in range(total_size)])
+
+    def reduce_sum_step_indexes(self, step, array):
+       return sum([array[i] for i in self.step_indexes(step)])
+
     def initialize(self):
         grid_array = [(self.sim.grid.min(d), self.sim.grid.max(d)) for d in range(self.sim.ndims())]
         Call_Void(self.sim, "pairs->initDomain", [param for delim in grid_array for param in delim])
@@ -96,7 +106,20 @@ class BlockForest:
         return 1
 
     def step_indexes(self, step):
-        return [step]
+        yield from For(self.sim, 0, self.nranks)
+
+    def first_step_index(self, step):
+        return 0
+
+    def reduce_sum_all_steps(self, array):
+        return self.reduce_sum_step_indexes(0, array)
+
+    def reduce_sum_step_indexes(self, step, array):
+        nsend_sum = self.sim.add_temp_var(0)
+        for i in For(self.sim, 0, self.nranks):
+            Assign(self.sim, nsend_sum, nsend_sum + array[i])
+
+        return nsend_sum
 
     def initialize(self):
         grid_array = [(self.sim.grid.min(d), self.sim.grid.max(d)) for d in range(self.sim.ndims())]