diff --git a/src/cell.c b/src/cell.c
index 7ced653d2e1da294a598ea86045e6e3bb5577f63..4502f5d265dc68540e16ed0e51e681cf5733f842 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1893,6 +1893,13 @@ void cell_set_super(struct cell *c, struct cell *super) {
       if (c->progeny[k] != NULL) cell_set_super(c->progeny[k], super);
 }
 
+void cell_set_super_mapper(void *map_data, int num_elements, void *extra_data) {
+  for (int ind = 0; ind < num_elements; ind++) {
+    struct cell *c = &((struct cell *)map_data)[ind];
+    cell_set_super(c, NULL);
+  }
+}
+
 /**
  * @brief Recursively drifts the #part in a cell hierarchy.
  *
diff --git a/src/cell.h b/src/cell.h
index 82035695d465a681539de86db7c3bf68cad928da..e97400623dbb7a66aee981d21883fe4d8f73406a 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -404,6 +404,7 @@ void cell_activate_subcell_tasks(struct cell *ci, struct cell *cj,
 void cell_activate_drift_part(struct cell *c, struct scheduler *s);
 void cell_activate_sorts(struct cell *c, int sid, struct scheduler *s);
 void cell_clear_drift_flags(struct cell *c, void *data);
+void cell_set_super_mapper(void *map_data, int num_elements, void *extra_data);
 
 /* Inlined functions (for speed). */
 
diff --git a/src/engine.c b/src/engine.c
index cd3d05c09d54f54a4ed042c3c04ce0a0485c4fd5..8c26099204a69a68af121d44776cf5d1ad60a813 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -253,6 +253,16 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
   }
 }
 
+void engine_make_hierarchical_tasks_mapper(void *map_data, int num_elements,
+                                           void *extra_data) {
+  struct engine *e = (struct engine *)extra_data;
+
+  for (int ind = 0; ind < num_elements; ind++) {
+    struct cell *c = &((struct cell *)map_data)[ind];
+    engine_make_hierarchical_tasks(e, c);
+  }
+}
+
 #ifdef WITH_MPI
 /**
  * Do the exchange of one type of particles with all the other nodes.
@@ -1672,6 +1682,83 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
 #endif
 }
 
+void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
+                                           void *extra_data) {
+
+  struct engine *e = ((struct engine **)extra_data)[0];
+  struct task **ghosts = ((struct task ***)extra_data)[1];
+
+  struct space *s = e->s;
+  struct scheduler *sched = &e->sched;
+  const int nodeID = e->nodeID;
+  const int periodic = s->periodic;
+  const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
+  const int cdim_ghost[3] = {s->cdim[0] / 4 + 1, s->cdim[1] / 4 + 1,
+                             s->cdim[2] / 4 + 1};
+  const double theta_crit_inv = e->gravity_properties->theta_crit_inv;
+  struct cell *cells = s->cells_top;
+  const int n_ghosts = cdim_ghost[0] * cdim_ghost[1] * cdim_ghost[2] * 2;
+
+  /* Loop through the elements, which are just byte offsets from NULL. */
+  for (int ind = 0; ind < num_elements; ind++) {
+
+    /* Get the cell index. */
+    const int cid = (size_t)(map_data) + ind;
+    const int i = cid / (cdim[1] * cdim[2]);
+    const int j = (cid / cdim[2]) % cdim[1];
+    const int k = cid % cdim[2];
+
+    /* Get the cell */
+    struct cell *ci = &cells[cid];
+
+    /* Skip cells without gravity particles */
+    if (ci->gcount == 0) continue;
+
+    /* Is that cell local ? */
+    if (ci->nodeID != nodeID) continue;
+
+    /* If the cells is local build a self-interaction */
+    scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci, NULL);
+
+    /* Deal with periodicity dependencies */
+    const int ghost_id = cell_getid(cdim_ghost, i / 4, j / 4, k / 4);
+    if (ghost_id > n_ghosts) error("Invalid ghost_id");
+    if (periodic) {
+      ci->grav_ghost[0] = ghosts[2 * ghost_id + 0];
+      ci->grav_ghost[1] = ghosts[2 * ghost_id + 1];
+    }
+
+    /* Loop over every other cell */
+    for (int ii = 0; ii < cdim[0]; ii++) {
+      for (int jj = 0; jj < cdim[1]; jj++) {
+        for (int kk = 0; kk < cdim[2]; kk++) {
+
+          /* Get the cell */
+          const int cjd = cell_getid(cdim, ii, jj, kk);
+          struct cell *cj = &cells[cjd];
+
+          /* Avoid duplicates */
+          if (cid <= cjd) continue;
+
+          /* Skip cells without gravity particles */
+          if (cj->gcount == 0) continue;
+
+          /* Is that neighbour local ? */
+          if (cj->nodeID != nodeID) continue;  // MATTHIEU
+
+          /* Are the cells to close for a MM interaction ? */
+          if (!gravity_multipole_accept(ci->multipole, cj->multipole,
+                                        theta_crit_inv, 1)) {
+
+            scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0,
+                              ci, cj);
+          }
+        }
+      }
+    }
+  }
+}
+
 /**
  * @brief Constructs the top-level tasks for the short-range gravity
  * interactions.
@@ -1686,13 +1773,9 @@ void engine_make_self_gravity_tasks(struct engine *e) {
 
   struct space *s = e->s;
   struct scheduler *sched = &e->sched;
-  const int nodeID = e->nodeID;
   const int periodic = s->periodic;
-  const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
   const int cdim_ghost[3] = {s->cdim[0] / 4 + 1, s->cdim[1] / 4 + 1,
                              s->cdim[2] / 4 + 1};
-  const double theta_crit_inv = e->gravity_properties->theta_crit_inv;
-  struct cell *cells = s->cells_top;
   struct task **ghosts = NULL;
   const int n_ghosts = cdim_ghost[0] * cdim_ghost[1] * cdim_ghost[2] * 2;
 
@@ -1720,64 +1803,12 @@ void engine_make_self_gravity_tasks(struct engine *e) {
     }
   }
 
-  /* Run through the higher level cells */
-  for (int i = 0; i < cdim[0]; i++) {
-    for (int j = 0; j < cdim[1]; j++) {
-      for (int k = 0; k < cdim[2]; k++) {
-
-        /* Get the cell */
-        const int cid = cell_getid(cdim, i, j, k);
-        struct cell *ci = &cells[cid];
-
-        /* Skip cells without gravity particles */
-        if (ci->gcount == 0) continue;
+  /* Cretae the multipole self and pair tasks. */
+  void *extra_data[2] = {e, ghosts};
+  threadpool_map(&e->threadpool, engine_make_self_gravity_tasks_mapper, NULL,
+                 s->nr_cells, 1, 0, extra_data);
 
-        /* Is that cell local ? */
-        if (ci->nodeID != nodeID) continue;
-
-        /* If the cells is local build a self-interaction */
-        scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci,
-                          NULL);
-
-        /* Deal with periodicity dependencies */
-        const int ghost_id = cell_getid(cdim_ghost, i / 4, j / 4, k / 4);
-        if (ghost_id > n_ghosts) error("Invalid ghost_id");
-        if (periodic) {
-          ci->grav_ghost[0] = ghosts[2 * ghost_id + 0];
-          ci->grav_ghost[1] = ghosts[2 * ghost_id + 1];
-        }
-
-        /* Loop over every other cell */
-        for (int ii = 0; ii < cdim[0]; ii++) {
-          for (int jj = 0; jj < cdim[1]; jj++) {
-            for (int kk = 0; kk < cdim[2]; kk++) {
-
-              /* Get the cell */
-              const int cjd = cell_getid(cdim, ii, jj, kk);
-              struct cell *cj = &cells[cjd];
-
-              /* Avoid duplicates */
-              if (cid <= cjd) continue;
-
-              /* Skip cells without gravity particles */
-              if (cj->gcount == 0) continue;
-
-              /* Is that neighbour local ? */
-              if (cj->nodeID != nodeID) continue;  // MATTHIEU
-
-              /* Are the cells to close for a MM interaction ? */
-              if (!gravity_multipole_accept(ci->multipole, cj->multipole,
-                                            theta_crit_inv, 1)) {
-
-                scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0,
-                                  0, ci, cj);
-              }
-            }
-          }
-        }
-      }
-    }
-  }
+  /* Clean up. */
   if (periodic) free(ghosts);
 }
 
@@ -1814,9 +1845,15 @@ void engine_make_external_gravity_tasks(struct engine *e) {
  * Additional loop over neighbours can later be added by simply duplicating
  * all the tasks created by this function.
  *
- * @param e The #engine.
+ * @param map_data Offset of first two indices disguised as a pointer.
+ * @param num_elements Number of cells to traverse.
+ * @param extra_data The #engine.
  */
-void engine_make_hydroloop_tasks(struct engine *e) {
+void engine_make_hydroloop_tasks_mapper(void *map_data, int num_elements,
+                                        void *extra_data) {
+
+  /* Extract the engine pointer. */
+  struct engine *e = (struct engine *)extra_data;
 
   struct space *s = e->s;
   struct scheduler *sched = &e->sched;
@@ -1824,53 +1861,53 @@ void engine_make_hydroloop_tasks(struct engine *e) {
   const int *cdim = s->cdim;
   struct cell *cells = s->cells_top;
 
-  /* Run through the highest level of cells and add pairs. */
-  for (int i = 0; i < cdim[0]; i++) {
-    for (int j = 0; j < cdim[1]; j++) {
-      for (int k = 0; k < cdim[2]; k++) {
-
-        /* Get the cell */
-        const int cid = cell_getid(cdim, i, j, k);
-        struct cell *ci = &cells[cid];
-
-        /* Skip cells without hydro particles */
-        if (ci->count == 0) continue;
-
-        /* If the cells is local build a self-interaction */
-        if (ci->nodeID == nodeID)
-          scheduler_addtask(sched, task_type_self, task_subtype_density, 0, 0,
-                            ci, NULL);
-
-        /* Now loop over all the neighbours of this cell */
-        for (int ii = -1; ii < 2; ii++) {
-          int iii = i + ii;
-          if (!s->periodic && (iii < 0 || iii >= cdim[0])) continue;
-          iii = (iii + cdim[0]) % cdim[0];
-          for (int jj = -1; jj < 2; jj++) {
-            int jjj = j + jj;
-            if (!s->periodic && (jjj < 0 || jjj >= cdim[1])) continue;
-            jjj = (jjj + cdim[1]) % cdim[1];
-            for (int kk = -1; kk < 2; kk++) {
-              int kkk = k + kk;
-              if (!s->periodic && (kkk < 0 || kkk >= cdim[2])) continue;
-              kkk = (kkk + cdim[2]) % cdim[2];
-
-              /* Get the neighbouring cell */
-              const int cjd = cell_getid(cdim, iii, jjj, kkk);
-              struct cell *cj = &cells[cjd];
-
-              /* Is that neighbour local and does it have particles ? */
-              if (cid >= cjd || cj->count == 0 ||
-                  (ci->nodeID != nodeID && cj->nodeID != nodeID))
-                continue;
-
-              /* Construct the pair task */
-              const int sid =
-                  sortlistID[(kk + 1) + 3 * ((jj + 1) + 3 * (ii + 1))];
-              scheduler_addtask(sched, task_type_pair, task_subtype_density,
-                                sid, 0, ci, cj);
-            }
-          }
+  /* Loop through the elements, which are just byte offsets from NULL. */
+  for (int ind = 0; ind < num_elements; ind++) {
+
+    /* Get the cell index. */
+    const int cid = (size_t)(map_data) + ind;
+    const int i = cid / (cdim[1] * cdim[2]);
+    const int j = (cid / cdim[2]) % cdim[1];
+    const int k = cid % cdim[2];
+
+    /* Get the cell */
+    struct cell *ci = &cells[cid];
+
+    /* Skip cells without hydro particles */
+    if (ci->count == 0) continue;
+
+    /* If the cells is local build a self-interaction */
+    if (ci->nodeID == nodeID)
+      scheduler_addtask(sched, task_type_self, task_subtype_density, 0, 0, ci,
+                        NULL);
+
+    /* Now loop over all the neighbours of this cell */
+    for (int ii = -1; ii < 2; ii++) {
+      int iii = i + ii;
+      if (!s->periodic && (iii < 0 || iii >= cdim[0])) continue;
+      iii = (iii + cdim[0]) % cdim[0];
+      for (int jj = -1; jj < 2; jj++) {
+        int jjj = j + jj;
+        if (!s->periodic && (jjj < 0 || jjj >= cdim[1])) continue;
+        jjj = (jjj + cdim[1]) % cdim[1];
+        for (int kk = -1; kk < 2; kk++) {
+          int kkk = k + kk;
+          if (!s->periodic && (kkk < 0 || kkk >= cdim[2])) continue;
+          kkk = (kkk + cdim[2]) % cdim[2];
+
+          /* Get the neighbouring cell */
+          const int cjd = cell_getid(cdim, iii, jjj, kkk);
+          struct cell *cj = &cells[cjd];
+
+          /* Is that neighbour local and does it have particles ? */
+          if (cid >= cjd || cj->count == 0 ||
+              (ci->nodeID != nodeID && cj->nodeID != nodeID))
+            continue;
+
+          /* Construct the pair task */
+          const int sid = sortlistID[(kk + 1) + 3 * ((jj + 1) + 3 * (ii + 1))];
+          scheduler_addtask(sched, task_type_pair, task_subtype_density, sid, 0,
+                            ci, cj);
         }
       }
     }
@@ -1883,17 +1920,16 @@ void engine_make_hydroloop_tasks(struct engine *e) {
  * For each hydrodynamic and gravity task, construct the links with
  * the corresponding cell.  Similarly, construct the dependencies for
  * all the sorting tasks.
- *
- * @param e The #engine.
  */
-void engine_count_and_link_tasks(struct engine *e) {
+void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
+                                        void *extra_data) {
 
+  struct engine *e = (struct engine *)extra_data;
   struct scheduler *const sched = &e->sched;
-  const int nr_tasks = sched->nr_tasks;
 
-  for (int ind = 0; ind < nr_tasks; ind++) {
+  for (int ind = 0; ind < num_elements; ind++) {
+    struct task *const t = &((struct task *)map_data)[ind];
 
-    struct task *const t = &sched->tasks[ind];
     struct cell *const ci = t->ci;
     struct cell *const cj = t->cj;
 
@@ -2149,18 +2185,17 @@ static inline void engine_make_hydro_loops_dependencies(struct scheduler *sched,
  * corresponding to the second hydro loop over neighbours.
  * With all the relevant tasks for a given cell available, we construct
  * all the dependencies for that cell.
- *
- * @param e The #engine.
  */
-void engine_make_extra_hydroloop_tasks(struct engine *e) {
+void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
+                                              void *extra_data) {
 
+  struct engine *e = (struct engine *)extra_data;
   struct scheduler *sched = &e->sched;
-  const int nr_tasks = sched->nr_tasks;
   const int nodeID = e->nodeID;
   const int with_cooling = (e->policy & engine_policy_cooling);
 
-  for (int ind = 0; ind < nr_tasks; ind++) {
-    struct task *t = &sched->tasks[ind];
+  for (int ind = 0; ind < num_elements; ind++) {
+    struct task *t = &((struct task *)map_data)[ind];
 
     /* Sort tasks depend on the drift of the cell. */
     if (t->type == task_type_sort && t->ci->nodeID == engine_rank) {
@@ -2175,7 +2210,7 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
       scheduler_addunlock(sched, t->ci->super->sorts, t);
 
 #ifdef EXTRA_HYDRO_LOOP
-      /* Start by constructing the task for the second  and third hydro loop */
+      /* Start by constructing the task for the second  and third hydro loop. */
       struct task *t2 = scheduler_addtask(
           sched, task_type_self, task_subtype_gradient, 0, 0, t->ci, NULL);
       struct task *t3 = scheduler_addtask(
@@ -2424,7 +2459,10 @@ void engine_maketasks(struct engine *e) {
   scheduler_reset(sched, s->tot_cells * engine_maxtaskspercell);
 
   /* Construct the firt hydro loop over neighbours */
-  if (e->policy & engine_policy_hydro) engine_make_hydroloop_tasks(e);
+  if (e->policy & engine_policy_hydro) {
+    threadpool_map(&e->threadpool, engine_make_hydroloop_tasks_mapper, NULL,
+                   s->nr_cells, 1, 0, e);
+  }
 
   /* Add the self gravity tasks. */
   if (e->policy & engine_policy_self_gravity) engine_make_self_gravity_tasks(e);
@@ -2475,20 +2513,23 @@ void engine_maketasks(struct engine *e) {
   /* Count the number of tasks associated with each cell and
      store the density tasks in each cell, and make each sort
      depend on the sorts of its progeny. */
-  engine_count_and_link_tasks(e);
+  threadpool_map(&e->threadpool, engine_count_and_link_tasks_mapper,
+                 sched->tasks, sched->nr_tasks, sizeof(struct task), 0, e);
 
   /* Now that the self/pair tasks are at the right level, set the super
    * pointers. */
-  for (int k = 0; k < nr_cells; k++) cell_set_super(&cells[k], NULL);
+  threadpool_map(&e->threadpool, cell_set_super_mapper, cells, nr_cells,
+                 sizeof(struct cell), 0, NULL);
 
   /* Append hierarchical tasks to each cell. */
-  for (int k = 0; k < nr_cells; k++)
-    engine_make_hierarchical_tasks(e, &cells[k]);
+  threadpool_map(&e->threadpool, engine_make_hierarchical_tasks_mapper, cells,
+                 nr_cells, sizeof(struct cell), 0, e);
 
   /* Run through the tasks and make force tasks for each density task.
      Each force task depends on the cell ghosts and unlocks the kick task
      of its super-cell. */
-  if (e->policy & engine_policy_hydro) engine_make_extra_hydroloop_tasks(e);
+  threadpool_map(&e->threadpool, engine_make_extra_hydroloop_tasks_mapper,
+                 sched->tasks, sched->nr_tasks, sizeof(struct task), 0, e);
 
   /* Add the dependencies for the gravity stuff */
   if (e->policy & (engine_policy_self_gravity | engine_policy_external_gravity))
@@ -2641,8 +2682,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
             scheduler_activate(s, l->t);
 
             /* Drift the cell which will be sent at the level at which it is
-               sent,
-               i.e. drift the cell specified in the send task (l->t) itself. */
+               sent, i.e. drift the cell specified in the send task (l->t)
+               itself. */
             cell_activate_drift_part(l->t->ci, s);
 
             if (cell_is_active(cj, e)) {
@@ -2699,8 +2740,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
             scheduler_activate(s, l->t);
 
             /* Drift the cell which will be sent at the level at which it is
-               sent,
-               i.e. drift the cell specified in the send task (l->t) itself. */
+               sent, i.e. drift the cell specified in the send task (l->t)
+               itself. */
             cell_activate_drift_part(l->t->ci, s);
 
             if (cell_is_active(ci, e)) {