Added ghost tasks to linke dependencies of the FFT mesh task

examples/UniformDMBox/makeIC.py

@@ -26,7 +26,7 @@ from numpy import *
 # with a density of 1
 
 # Parameters
-periodic= 0           # 1 For periodic box
+periodic= 1           # 1 For periodic box
 boxSize = 1.
 rho = 1.
 L = int(sys.argv[1])  # Number of particles along one axis

examples/UniformDMBox/uniformBox.yml

@@ -14,7 +14,7 @@ TimeIntegration:
   dt_max:     1e-3  # The maximal time-step size of the simulation (in internal units).
 
 Scheduler:
-  max_top_level_cells: 16
+  max_top_level_cells: 8
   cell_split_size:     50
   
 # Parameters governing the snapshots

src/cell.c

@@ -1467,6 +1467,28 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
     }
   }
 
+  /* for (struct link *l = c->grav; l != NULL; l = l->next) { */
+
+  /*   struct task *t = l->t; */
+  /*   struct cell *ci = t->ci; */
+  /*   struct cell *cj = t->cj; */
+  /*   scheduler_activate(s, t); */
+
+  /*   if (t->type == task_type_pair) { */
+  /*     scheduler_activate(s, ci->drift_gpart); */
+  /*     scheduler_activate(s, cj->drift_gpart); */
+  /*     scheduler_activate(s, ci->init_grav); */
+  /*     scheduler_activate(s, cj->init_grav); */
+  /*     scheduler_activate(s, ci->grav_ghost[1]); */
+  /*     scheduler_activate(s, cj->grav_ghost[1]); */
+  /*   } */
+
+  /*   if (t->type == task_type_self) { */
+  /*     scheduler_activate(s, ci->drift_gpart); */
+  /*     scheduler_activate(s, ci->init_grav); */
+  /*   } */
+  /* } */
+
   /* Unskip all the other task types. */
   for (struct link *l = c->gradient; l != NULL; l = l->next)
     scheduler_activate(s, l->t);
@@ -1482,6 +1504,8 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
   if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
   if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
   if (c->timestep != NULL) scheduler_activate(s, c->timestep);
+  if (c->grav_ghost[0] != NULL) scheduler_activate(s, c->grav_ghost[0]);
+  if (c->grav_ghost[1] != NULL) scheduler_activate(s, c->grav_ghost[1]);
   if (c->grav_down != NULL) scheduler_activate(s, c->grav_down);
   if (c->grav_long_range != NULL) scheduler_activate(s, c->grav_long_range);
   if (c->cooling != NULL) scheduler_activate(s, c->cooling);

src/cell.h

@@ -172,6 +172,9 @@ struct cell {
   /*! The task to compute time-steps */
   struct task *timestep;
 
+  /*! Task linking the FFT mesh to the rest of gravity tasks */
+  struct task *grav_ghost[2];
+
   /*! Task computing long range non-periodic gravity interactions */
   struct task *grav_long_range;
 

src/engine.c

@@ -133,6 +133,7 @@ void engine_addlink(struct engine *e, struct link **l, struct task *t) {
 void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
 
   struct scheduler *s = &e->sched;
+  const int periodic = e->s->periodic;
   const int is_hydro = (e->policy & engine_policy_hydro);
   const int is_self_gravity = (e->policy & engine_policy_self_gravity);
   const int is_with_cooling = (e->policy & engine_policy_cooling);
@@ -172,6 +173,7 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
         c->grav_down = scheduler_addtask(s, task_type_grav_down,
                                          task_subtype_none, 0, 0, c, NULL);
 
+        if (periodic) scheduler_addunlock(s, c->init_grav, c->grav_ghost[0]);
         scheduler_addunlock(s, c->init_grav, c->grav_long_range);
         scheduler_addunlock(s, c->grav_long_range, c->grav_down);
         scheduler_addunlock(s, c->grav_down, c->kick2);
@@ -1645,51 +1647,108 @@ 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;
-  const int nr_cells = s->nr_cells;
+  struct task **ghosts = NULL;
+  const int n_ghosts = cdim_ghost[0] * cdim_ghost[1] * cdim_ghost[2] * 2;
+
+  /* Create the top-level task if periodic */
+  if (periodic) {
+
+    /* Create the FFT task for this MPI rank */
+    s->grav_top_level = scheduler_addtask(sched, task_type_grav_top_level,
+                                          task_type_none, 0, 0, NULL, NULL);
+
+    /* Create a grid of ghosts to deal with the dependencies */
+    if ((ghosts = malloc(n_ghosts * sizeof(struct task *))) == 0)
+      error("Error allocating memory for gravity fft ghosts");
+
+    /* Make the ghosts implicit and add the dependencies */
+    for (int n = 0; n < n_ghosts / 2; ++n) {
+      ghosts[2 * n + 0] = scheduler_addtask(sched, task_type_grav_ghost,
+                                            task_type_none, 0, 0, NULL, NULL);
+      ghosts[2 * n + 1] = scheduler_addtask(sched, task_type_grav_ghost,
+                                            task_type_none, 0, 0, NULL, NULL);
+      ghosts[2 * n + 0]->implicit = 1;
+      ghosts[2 * n + 1]->implicit = 1;
+      scheduler_addunlock(sched, ghosts[2 * n + 0], s->grav_top_level);
+      scheduler_addunlock(sched, s->grav_top_level, ghosts[2 * n + 1]);
+    }
+  }
 
-  /* Create the FFT task for this MPI rank */
-  struct task *fft = scheduler_addtask(sched, task_type_grav_top_level,
-                                       task_type_none, 0, 0, NULL, NULL);
+  /* 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++) {
 
-  for (int cid = 0; cid < nr_cells; ++cid) {
+        /* Get the cell */
+        const int cid = cell_getid(cdim, i, j, k);
+        struct cell *ci = &cells[cid];
 
-    struct cell *ci = &cells[cid];
+        /* Skip cells without gravity particles */
+        if (ci->gcount == 0) continue;
 
-    /* Skip cells without gravity particles */
-    if (ci->gcount == 0) continue;
+        /* Is that cell local ? */
+        if (ci->nodeID != nodeID) continue;
 
-    /* Is that cell local ? */
-    if (ci->nodeID != nodeID) continue;
+        /* If the cells is local build a self-interaction */
+        struct task *self = 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) {
+          scheduler_addunlock(sched, ghosts[2 * ghost_id + 1], self);
+          ci->grav_ghost[0] = ghosts[2 * ghost_id + 0];
+          ci->grav_ghost[1] = ghosts[2 * ghost_id + 1];
+        }
 
-    /* If the cells is local build a self-interaction */
-    struct task *self = scheduler_addtask(sched, task_type_self,
-                                          task_subtype_grav, 0, 0, ci, NULL);
+        /* 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++) {
 
-    scheduler_addunlock(sched, fft, self);
+              /* Get the cell */
+              const int cjd = cell_getid(cdim, ii, jj, kk);
+              const int ghost_jd =
+                  cell_getid(cdim_ghost, ii / 4, jj / 4, kk / 4);
+              struct cell *cj = &cells[cjd];
 
-    /* Loop over every other cell */
-    for (int cjd = cid + 1; cjd < nr_cells; ++cjd) {
+              /* Avoid duplicates */
+              if (cid <= cjd) continue;
 
-      struct cell *cj = &cells[cjd];
+              /* Skip cells without gravity particles */
+              if (cj->gcount == 0) continue;
 
-      /* Skip cells without gravity particles */
-      if (cj->gcount == 0) continue;
+              /* Is that neighbour local ? */
+              if (cj->nodeID != nodeID) continue;  // MATTHIEU
 
-      /* 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)) {
 
-      /* 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);
+                struct task *pair = scheduler_addtask(
+                    sched, task_type_pair, task_subtype_grav, 0, 0, ci, cj);
 
-        // scheduler_addunlock(sched, fft, pair);
+                /* Deal with periodicity dependencies */
+                if (periodic) {
+                  scheduler_addunlock(sched, ghosts[2 * ghost_id + 1], pair);
+                  if (ghost_id != ghost_jd)
+                    scheduler_addunlock(sched, ghosts[2 * ghost_jd + 1], pair);
+                }
+              }
+            }
+          }
+        }
       }
     }
   }
+  if (periodic) free(ghosts);
 }
 
 void engine_make_external_gravity_tasks(struct engine *e) {
@@ -2618,7 +2677,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
     }
 
     /* Periodic gravity ? */
-    else if (t->type == task_type_grav_top_level) {
+    else if (t->type == task_type_grav_top_level ||
+             t->type == task_type_grav_ghost) {
       scheduler_activate(s, t);
     }
 
@@ -3053,6 +3113,7 @@ void engine_skip_force_and_kick(struct engine *e) {
         t->type == task_type_timestep || t->subtype == task_subtype_force ||
         t->subtype == task_subtype_grav ||
         t->type == task_type_grav_long_range ||
+        t->type == task_type_grav_ghost ||
         t->type == task_type_grav_top_level || t->type == task_type_grav_down ||
         t->type == task_type_cooling || t->type == task_type_sourceterms)
       t->skip = 1;
@@ -3304,8 +3365,8 @@ void engine_step(struct engine *e) {
 
     if (e->policy & engine_policy_reconstruct_mpoles)
       engine_reconstruct_multipoles(e);
-    else
-      engine_drift_top_multipoles(e);
+    // else
+    //  engine_drift_top_multipoles(e);
   }
 
   /* Print the number of active tasks ? */
@@ -3417,9 +3478,14 @@ int engine_is_done(struct engine *e) {
 void engine_unskip(struct engine *e) {
 
   const ticks tic = getticks();
+
+  /* Activate all the regular tasks */
   threadpool_map(&e->threadpool, runner_do_unskip_mapper, e->s->cells_top,
                  e->s->nr_cells, sizeof(struct cell), 1, e);
 
+  /* And the top level gravity FFT one */
+  if (e->s->periodic) scheduler_activate(&e->sched, e->s->grav_top_level);
+
   if (e->verbose)
     message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
             clocks_getunit());

src/runner_doiact_fft.c

@@ -169,10 +169,10 @@ void runner_do_grav_fft(struct runner* r, int timer) {
   struct cell* cells = s->cells_top;
 
   /* Make sure everything has been drifted to the current point */
-  for (int i = 0; i < nr_cells; ++i) {
+  for (int i = 0; i < nr_cells; ++i)
     if (cells[i].ti_old_multipole != ti_current)
       cell_drift_multipole(&cells[i], e);
-  }
+  // error("Top-level multipole %d not drifted", i);
 
   /* Allocates some memory for the density mesh */
   double* restrict rho = fftw_alloc_real(N * N * N);

src/runner_doiact_grav.h

@@ -429,6 +429,11 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
 
   /* Sanity check */
   if (ci == cj) error("Pair interaction between a cell and itself.");
+
+  if (cell_is_active(ci, e) && ci->ti_old_multipole != e->ti_current)
+    error("ci->multipole not drifted.");
+  if (cell_is_active(cj, e) && cj->ti_old_multipole != e->ti_current)
+    error("cj->multipole not drifted.");
 #endif
 
 #if ICHECK > 0

src/scheduler.c

@@ -137,7 +137,8 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
     redo = 0;
 
     /* Non-splittable task? */
-    if ((t->ci == NULL && t->type != task_type_grav_top_level) ||
+    if ((t->ci == NULL && t->type != task_type_grav_top_level &&
+         t->type != task_type_grav_ghost) ||
         ((t->type == task_type_pair && t->cj == NULL)) ||
         ((t->type == task_type_kick1) && t->ci->nodeID != s->nodeID) ||
         ((t->type == task_type_kick2) && t->ci->nodeID != s->nodeID) ||
@@ -791,7 +792,8 @@ void scheduler_set_unlocks(struct scheduler *s) {
     /* Check that we are not overflowing */
     if (counts[s->unlock_ind[k]] < 0)
       error("Task unlocking more than %d other tasks!",
-            (1 << (sizeof(short int) - 1)) - 1);
+            (1 << (8 * sizeof(short int) - 1)) - 1);
+
 #endif
   }
 
@@ -1148,7 +1150,9 @@ void scheduler_start(struct scheduler *s) {
       if (t->type == task_type_send || t->type == task_type_recv) continue;
 
       /* Don't check the FFT task */
-      if (t->type == task_type_grav_top_level) continue;
+      if (t->type == task_type_grav_top_level ||
+          t->type == task_type_grav_ghost)
+        continue;
 
       if (ci == NULL && cj == NULL) {
 
@@ -1195,6 +1199,8 @@ void scheduler_start(struct scheduler *s) {
   }
 #endif
 
+  scheduler_print_tasks(s, "tasks.dat");
+
   /* Loop over the tasks and enqueue whoever is ready. */
   if (s->active_count > 1000) {
     threadpool_map(s->threadpool, scheduler_enqueue_mapper, s->tid_active,

src/space.c

@@ -223,6 +223,8 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
     c->drift_gpart = NULL;
     c->cooling = NULL;
     c->sourceterms = NULL;
+    c->grav_ghost[0] = NULL;
+    c->grav_ghost[1] = NULL;
     c->grav_long_range = NULL;
     c->grav_down = NULL;
     c->super = c;
@@ -2590,6 +2592,8 @@ void space_init(struct space *s, const struct swift_params *params,
   s->sparts = sparts;
   s->nr_queues = 1; /* Temporary value until engine construction */
 
+  s->periodic = 0;
+
   /* Are we replicating the space ? */
   if (replicate < 1)
     error("Value of 'InitialConditions:replicate' (%d) is too small",

src/space.h

@@ -130,6 +130,9 @@ struct space {
   /*! The s-particle data (cells have pointers to this). */
   struct spart *sparts;
 
+  /*! The top-level FFT task */
+  struct task *grav_top_level;
+
   /*! General-purpose lock for this space. */
   swift_lock_type lock;
 

src/task.c

@@ -54,8 +54,8 @@ const char *taskID_names[task_type_count] = {
     "drift_part", "drift_gpart",    "kick1",
     "kick2",      "timestep",       "send",
     "recv",       "grav_top_level", "grav_long_range",
-    "grav_mm",    "grav_down",      "cooling",
-    "sourceterms"};
+    "grav_ghost", "grav_mm",        "grav_down",
+    "cooling",    "sourceterms"};
 
 /* Sub-task type names. */
 const char *subtaskID_names[task_subtype_count] = {

src/task.h

@@ -56,6 +56,7 @@ enum task_types {
   task_type_recv,
   task_type_grav_top_level,
   task_type_grav_long_range,
+  task_type_grav_ghost,
   task_type_grav_mm,
   task_type_grav_down,
   task_type_cooling,