The top-level gravity task can now be run without being associated with a cell.

examples/EAGLE_12/eagle_12.yml

@@ -31,8 +31,6 @@ Gravity:
   eta:                   0.025    # Constant dimensionless multiplier for time integration.
   theta:                 0.7      # Opening angle (Multipole acceptance criterion)
   epsilon:               0.0001   # Softening length (in internal units).
-  a_smooth:              1000.
-  r_cut:                 4.
   
 # Parameters for the hydrodynamics scheme
 SPH:

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: 8
+  max_top_level_cells: 16
   cell_split_size:     50
   
 # Parameters governing the snapshots
@@ -35,4 +35,4 @@ Statistics:
 
 # Parameters related to the initial conditions
 InitialConditions:
-  file_name:  ./uniformDMBox_100.hdf5     # The file to read
+  file_name:  ./uniformDMBox_50.hdf5     # The file to read

src/cell.c

@@ -1453,7 +1453,6 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
   if (c->timestep != NULL) scheduler_activate(s, c->timestep);
   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->grav_top_level != NULL) scheduler_activate(s, c->grav_top_level);
   if (c->cooling != NULL) scheduler_activate(s, c->cooling);
   if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms);
 

src/cell.h

@@ -169,10 +169,7 @@ struct cell {
   /*! The task to compute time-steps */
   struct task *timestep;
 
-  /*! Task constructing the multipole from the particles */
-  struct task *grav_top_level;
-
-  /*! Task constructing the multipole from the particles */
+  /*! Task computing long range non-periodic gravity interactions */
   struct task *grav_long_range;
 
   /*! Task propagating the multipole to the particles */

src/engine.c

@@ -168,18 +168,12 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
         c->grav_long_range = scheduler_addtask(
             s, task_type_grav_long_range, task_subtype_none, 0, 0, c, NULL);
 
-        /* Gravity top-level periodic calculation */
-        c->grav_top_level = scheduler_addtask(s, task_type_grav_top_level,
-                                              task_subtype_none, 0, 0, c, NULL);
-
         /* Gravity recursive down-pass */
         c->grav_down = scheduler_addtask(s, task_type_grav_down,
                                          task_subtype_none, 0, 0, c, NULL);
 
         scheduler_addunlock(s, c->init_grav, c->grav_long_range);
-        scheduler_addunlock(s, c->init_grav, c->grav_top_level);
         scheduler_addunlock(s, c->grav_long_range, c->grav_down);
-        scheduler_addunlock(s, c->grav_top_level, c->grav_down);
         scheduler_addunlock(s, c->grav_down, c->kick2);
       }
 
@@ -1655,6 +1649,10 @@ void engine_make_self_gravity_tasks(struct engine *e) {
   struct cell *cells = s->cells_top;
   const int nr_cells = s->nr_cells;
 
+  /* 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);
+
   for (int cid = 0; cid < nr_cells; ++cid) {
 
     struct cell *ci = &cells[cid];
@@ -1666,7 +1664,10 @@ void engine_make_self_gravity_tasks(struct engine *e) {
     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);
+    struct task *self = scheduler_addtask(sched, task_type_self,
+                                          task_subtype_grav, 0, 0, ci, NULL);
+
+    scheduler_addunlock(sched, fft, self);
 
     /* Loop over every other cell */
     for (int cjd = cid + 1; cjd < nr_cells; ++cjd) {
@@ -1681,9 +1682,12 @@ void engine_make_self_gravity_tasks(struct engine *e) {
 
       /* Are the cells to close for a MM interaction ? */
       if (!gravity_multipole_accept(ci->multipole, cj->multipole,
-                                    theta_crit_inv, 1))
+                                    theta_crit_inv, 1)) {
         scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, ci,
                           cj);
+
+        // scheduler_addunlock(sched, fft, pair);
+      }
     }
   }
 }
@@ -2605,11 +2609,15 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
 
     /* Gravity ? */
     else if (t->type == task_type_grav_down ||
-             t->type == task_type_grav_long_range ||
-             t->type == task_type_grav_top_level) {
+             t->type == task_type_grav_long_range) {
       if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
     }
 
+    /* Periodic gravity ? */
+    else if (t->type == task_type_grav_top_level) {
+      scheduler_activate(s, t);
+    }
+
     /* Time-step? */
     else if (t->type == task_type_timestep) {
       t->ci->updated = 0;

src/gravity_properties.c

@@ -69,11 +69,9 @@ void gravity_props_print(const struct gravity_props *p) {
   message("Self-gravity softening:    epsilon=%.4f (Plummer equivalent: %.4f)",
           p->epsilon, p->epsilon / 3.);
 
-  if (p->a_smooth != gravity_props_default_a_smooth)
-    message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth);
+  message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth);
 
-  if (p->r_cut != gravity_props_default_r_cut)
-    message("Self-gravity MM cut-off: r_cut=%f", p->r_cut);
+  message("Self-gravity MM cut-off: r_cut=%f", p->r_cut);
 }
 
 #if defined(HAVE_HDF5)

src/runner.c

@@ -1728,10 +1728,12 @@ void *runner_main(void *data) {
 #ifdef SWIFT_DEBUG_CHECKS
       t->ti_run = e->ti_current;
 #ifndef WITH_MPI
-      if (ci == NULL && cj == NULL) {
+      if (t->type == task_type_grav_top_level) {
+        if (ci != NULL || cj != NULL)
+          error("Top-level gravity task associated with a cell");
+      } else if (ci == NULL && cj == NULL) {
 
         error("Task not associated with cells!");
-
       } else if (cj == NULL) { /* self */
 
         if (!cell_is_active(ci, e) && t->type != task_type_sort &&

src/runner_doiact_fft.c

@@ -47,7 +47,7 @@
  */
 __attribute__((always_inline)) INLINE static int row_major_id(int i, int j,
                                                               int k, int N) {
-  return (i * N * N + j * N + k);
+  return ((i % N) * N * N + (j % N) * N + (k % N));
 }
 
 /**
@@ -76,6 +76,12 @@ __attribute__((always_inline)) INLINE static void multipole_to_mesh_CIC(
   const double dz = fac * m->CoM[2] - k;
   const double tz = 1. - dz;
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (i < 0 || i >= N) error("Invalid multipole position in x");
+  if (j < 0 || j >= N) error("Invalid multipole position in y");
+  if (k < 0 || k >= N) error("Invalid multipole position in z");
+#endif
+
   /* CIC ! */
   rho[row_major_id(i + 0, j + 0, k + 0, N)] += m->m_pole.M_000 * tx * ty * tz;
   rho[row_major_id(i + 0, j + 0, k + 1, N)] += m->m_pole.M_000 * tx * ty * dz;
@@ -114,6 +120,12 @@ __attribute__((always_inline)) INLINE static void mesh_to_multipole_CIC(
   const double dz = fac * m->CoM[2] - k;
   const double tz = 1. - dz;
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (i < 0 || i >= N) error("Invalid multipole position in x");
+  if (j < 0 || j >= N) error("Invalid multipole position in y");
+  if (k < 0 || k >= N) error("Invalid multipole position in z");
+#endif
+
   /* CIC ! */
   m->pot.F_000 += pot[row_major_id(i + 0, j + 0, k + 0, N)] * tx * ty * tz;
   m->pot.F_000 += pot[row_major_id(i + 0, j + 0, k + 1, N)] * tx * ty * dz;
@@ -149,7 +161,6 @@ void runner_do_grav_fft(struct runner* r, int timer) {
   /* Some useful constants */
   const int N = cdim[0];
   const int N_half = N / 2;
-  const int Ncells = N * N * N;
   const double cell_fac = N / box_size;
 
   /* Recover the list of top-level multipoles */
@@ -164,7 +175,7 @@ void runner_do_grav_fft(struct runner* r, int timer) {
   }
 
   /* Allocates some memory for the density mesh */
-  double* restrict rho = fftw_alloc_real(Ncells);
+  double* restrict rho = fftw_alloc_real(N * N * N);
   if (rho == NULL) error("Error allocating memory for density mesh");
 
   /* Allocates some memory for the mesh in Fourier space */
@@ -179,7 +190,7 @@ void runner_do_grav_fft(struct runner* r, int timer) {
       N, N, N, frho, rho, FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
 
   /* Do a CIC mesh assignment of the multipoles */
-  bzero(rho, Ncells * sizeof(double));
+  bzero(rho, N * N * N * sizeof(double));
   for (int i = 0; i < nr_cells; ++i)
     multipole_to_mesh_CIC(&multipoles[i], rho, N, cell_fac);
 

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_pair && t->cj == NULL)) ||
+    if ((t->ci == NULL && t->type != task_type_grav_top_level) ||
+        ((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) ||
         ((t->type == task_type_drift) && t->ci->nodeID != s->nodeID) ||
@@ -1133,6 +1134,9 @@ void scheduler_start(struct scheduler *s) {
       /* Don't check MPI stuff */
       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 (ci == NULL && cj == NULL) {
 
         error("Task not associated with cells!");

src/space.c

@@ -26,6 +26,7 @@
 
 /* Some standard headers. */
 #include <float.h>
+#include <hdf5.h>
 #include <limits.h>
 #include <math.h>
 #include <stdlib.h>
@@ -220,7 +221,6 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
     c->drift = NULL;
     c->cooling = NULL;
     c->sourceterms = NULL;
-    c->grav_top_level = NULL;
     c->grav_long_range = NULL;
     c->grav_down = NULL;
     c->super = c;
@@ -2943,3 +2943,47 @@ void space_clean(struct space *s) {
   free(s->gparts);
   free(s->sparts);
 }
+
+void space_dump_multipoles(struct space *s) {
+
+  hid_t h_file =
+      H5Fcreate("multipoles.hdf5", H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+  hid_t h_grp =
+      H5Gcreate(h_file, "/Multipoles", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+
+  const hid_t h_space = H5Screate(H5S_SIMPLE);
+  const hid_t h_space2 = H5Screate(H5S_SIMPLE);
+  hsize_t shape[2] = {s->nr_cells, 3};
+  hsize_t shape2[1] = {s->nr_cells};
+  H5Sset_extent_simple(h_space, 2, shape, NULL);
+  H5Sset_extent_simple(h_space2, 1, shape2, NULL);
+
+  const hid_t h_data =
+      H5Dcreate(h_grp, "Coordinates", H5T_NATIVE_DOUBLE, h_space, H5P_DEFAULT,
+                H5P_DEFAULT, H5P_DEFAULT);
+
+  const hid_t h_data2 =
+      H5Dcreate(h_grp, "Potential", H5T_NATIVE_DOUBLE, h_space2, H5P_DEFAULT,
+                H5P_DEFAULT, H5P_DEFAULT);
+
+  double *temp = malloc(s->nr_cells * 3 * sizeof(double));
+  double *temp2 = malloc(s->nr_cells * 1 * sizeof(double));
+  for (int i = 0; i < s->nr_cells; ++i) {
+    temp[i * 3 + 0] = s->multipoles_top[i].CoM[0];
+    temp[i * 3 + 1] = s->multipoles_top[i].CoM[1];
+    temp[i * 3 + 2] = s->multipoles_top[i].CoM[2];
+    temp2[i] = s->multipoles_top[i].pot.F_000;
+  }
+
+  H5Dwrite(h_data, H5T_NATIVE_DOUBLE, h_space, H5S_ALL, H5P_DEFAULT, temp);
+  H5Dwrite(h_data2, H5T_NATIVE_DOUBLE, h_space2, H5S_ALL, H5P_DEFAULT, temp2);
+
+  free(temp);
+  free(temp2);
+  H5Dclose(h_data);
+  H5Dclose(h_data2);
+  H5Sclose(h_space);
+  H5Sclose(h_space2);
+  H5Gclose(h_grp);
+  H5Fclose(h_file);
+}

src/space.h

@@ -221,5 +221,6 @@ void space_check_timesteps(struct space *s);
 void space_replicate(struct space *s, int replicate, int verbose);
 void space_reset_task_counters(struct space *s);
 void space_clean(struct space *s);
+void space_dump_multipoles(struct space *s);
 
 #endif /* SWIFT_SPACE_H */

src/task.c

@@ -323,7 +323,6 @@ void task_unlock(struct task *t) {
       cell_munlocktree(ci);
       break;
 
-    case task_type_grav_top_level:
     case task_type_grav_long_range:
     case task_type_grav_mm:
       cell_munlocktree(ci);
@@ -442,7 +441,6 @@ int task_lock(struct task *t) {
       }
       break;
 
-    case task_type_grav_top_level:
     case task_type_grav_long_range:
     case task_type_grav_mm:
       /* Lock the m-poles */