Merge branch 'gravity_multi_dt' into 'master'

Split the drift into hydro and gravity drift tasks. Implement FFT task and dependencies

Changes:
 - Two drift tasks, one pert particle type.
 - FFT task with gravity ghost tasks for the dependencies.

It passes my usual tests but could you make sure it does not break anything ? That should be the last big change. The rest of gravity should be contained within the tasks themselves.

See merge request !340

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/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

@@ -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/Makefile.am

@@ -63,7 +63,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
 		 kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h runner_doiact_fft.h \
                  runner_doiact_nosort.h units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
-		 dimension.h equation_of_state.h part_type.h \
+		 dimension.h equation_of_state.h part_type.h periodic.h \
 		 gravity.h gravity_io.h \
 		 gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
 		 gravity/Default/gravity_debug.h gravity/Default/gravity_part.h  \

src/active.h

@@ -29,25 +29,48 @@
 #include "timeline.h"
 
 /**
- * @brief Check that a cell been drifted to the current time.
+ * @brief Check that the #part in a #cell have been drifted to the current time.
  *
  * @param c The #cell.
  * @param e The #engine containing information about the current time.
  * @return 1 if the #cell has been drifted to the current time, 0 otherwise.
  */
-__attribute__((always_inline)) INLINE static int cell_is_drifted(
+__attribute__((always_inline)) INLINE static int cell_are_part_drifted(
     const struct cell *c, const struct engine *e) {
 
 #ifdef SWIFT_DEBUG_CHECKS
-  if (c->ti_old > e->ti_current)
+  if (c->ti_old_part > e->ti_current)
     error(
         "Cell has been drifted too far forward in time! c->ti_old=%lld (t=%e) "
         "and e->ti_current=%lld (t=%e)",
-        c->ti_old, c->ti_old * e->timeBase, e->ti_current,
+        c->ti_old_part, c->ti_old_part * e->timeBase, e->ti_current,
         e->ti_current * e->timeBase);
 #endif
 
-  return (c->ti_old == e->ti_current);
+  return (c->ti_old_part == e->ti_current);
+}
+
+/**
+ * @brief Check that the #gpart in a #cell have been drifted to the current
+ * time.
+ *
+ * @param c The #cell.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #cell has been drifted to the current time, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int cell_are_gpart_drifted(
+    const struct cell *c, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (c->ti_old_gpart > e->ti_current)
+    error(
+        "Cell has been drifted too far forward in time! c->ti_old=%lld (t=%e) "
+        "and e->ti_current=%lld (t=%e)",
+        c->ti_old_gpart, c->ti_old_gpart * e->timeBase, e->ti_current,
+        e->ti_current * e->timeBase);
+#endif
+
+  return (c->ti_old_gpart == e->ti_current);
 }
 
 /* Are cells / particles active for regular tasks ? */

src/cell.h

@@ -74,7 +74,7 @@ struct pcell {
 
   /* Stats on this cell's particles. */
   double h_max;
-  integertime_t ti_end_min, ti_end_max, ti_beg_max, ti_old;
+  integertime_t ti_end_min, ti_end_max, ti_beg_max, ti_old_part, ti_old_gpart;
 
   /* Number of particles in this cell. */
   int count, gcount, scount;
@@ -157,8 +157,11 @@ struct cell {
   /*! The extra ghost task for complex hydro schemes */
   struct task *extra_ghost;
 
-  /*! The drift task */
-  struct task *drift;
+  /*! The drift task for parts */
+  struct task *drift_part;
+
+  /*! The drift task for gparts */
+  struct task *drift_gpart;
 
   /*! The first kick task */
   struct task *kick1;
@@ -169,10 +172,10 @@ struct cell {
   /*! The task to compute time-steps */
   struct task *timestep;
 
-  /*! Task constructing the multipole from the particles */
-  struct task *grav_top_level;
+  /*! Task linking the FFT mesh to the rest of gravity tasks */
+  struct task *grav_ghost[2];
 
-  /*! 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 */
@@ -233,24 +236,30 @@ struct cell {
   /*! Maximum beginning of (integer) time step in this cell. */
   integertime_t ti_beg_max;
 
-  /*! Last (integer) time the cell's particle was drifted forward in time. */
-  integertime_t ti_old;
-
   /*! Last (integer) time the cell's sort arrays were updated. */
   integertime_t ti_sort;
 
+  /*! Last (integer) time the cell's part were drifted forward in time. */
+  integertime_t ti_old_part;
+
+  /*! Last (integer) time the cell's gpart were drifted forward in time. */
+  integertime_t ti_old_gpart;
+
   /*! Last (integer) time the cell's multipole was drifted forward in time. */
   integertime_t ti_old_multipole;
 
   /*! Minimum dimension, i.e. smallest edge of this cell (min(width)). */
   float dmin;
 
-  /*! Maximum particle movement in this cell since last construction. */
-  float dx_max;
-
   /*! Maximum particle movement in this cell since the last sort. */
   float dx_max_sort;
 
+  /*! Maximum part movement in this cell since last construction. */
+  float dx_max_part;
+
+  /*! Maximum gpart movement in this cell since last construction. */
+  float dx_max_gpart;
+
   /*! Nr of #part in this cell. */
   int count;
 
@@ -357,13 +366,15 @@ void cell_clean_links(struct cell *c, void *data);
 void cell_make_multipoles(struct cell *c, integertime_t ti_current);
 void cell_check_multipole(struct cell *c, void *data);
 void cell_clean(struct cell *c);
-void cell_check_particle_drift_point(struct cell *c, void *data);
+void cell_check_part_drift_point(struct cell *c, void *data);
+void cell_check_gpart_drift_point(struct cell *c, void *data);
 void cell_check_multipole_drift_point(struct cell *c, void *data);
 void cell_reset_task_counters(struct cell *c);
 int cell_is_drift_needed(struct cell *c, const struct engine *e);
 int cell_unskip_tasks(struct cell *c, struct scheduler *s);
 void cell_set_super(struct cell *c, struct cell *super);
-void cell_drift_particles(struct cell *c, const struct engine *e);
+void cell_drift_part(struct cell *c, const struct engine *e);
+void cell_drift_gpart(struct cell *c, const struct engine *e);
 void cell_drift_multipole(struct cell *c, const struct engine *e);
 void cell_drift_all_multipoles(struct cell *c, const struct engine *e);
 void cell_check_timesteps(struct cell *c);

src/debug.c

@@ -259,8 +259,8 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
     message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
     result = 0;
   }
-  if (c->dx_max != dx_max) {
-    message("%d Inconsistent dx_max: %f != %f", *depth, c->dx_max, dx_max);
+  if (c->dx_max_part != dx_max) {
+    message("%d Inconsistent dx_max: %f != %f", *depth, c->dx_max_part, dx_max);
     message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
     result = 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/minmax.h

@@ -43,18 +43,4 @@
     _a > _b ? _a : _b;            \
   })
 
-/**
- * @brief Limits the value of x to be between a and b
- *
- * Only wraps once. If x > 2b, the returned value will be larger than b.
- * Similarly for x < -b.
- */
-#define box_wrap(x, a, b)                               \
-  ({                                                    \
-    const __typeof__(x) _x = (x);                       \
-    const __typeof__(a) _a = (a);                       \
-    const __typeof__(b) _b = (b);                       \
-    _x < _a ? (_x + _b) : ((_x > _b) ? (_x - _b) : _x); \
-  })
-
 #endif /* SWIFT_MINMAX_H */

src/multipole.h

@@ -28,7 +28,6 @@
 #include <string.h>
 
 /* Includes. */
-//#include "active.h"
 #include "align.h"
 #include "const.h"
 #include "error.h"
@@ -37,8 +36,8 @@
 #include "gravity_softened_derivatives.h"
 #include "inline.h"
 #include "kernel_gravity.h"
-#include "minmax.h"
 #include "part.h"
+#include "periodic.h"
 #include "vector_power.h"
 
 #define multipole_align 128

src/partition.c

@@ -524,7 +524,7 @@ static void repart_edge_metis(int partweights, int bothweights, int nodeID,
         t->type != task_type_sub_self && t->type != task_type_sub_self &&
         t->type != task_type_ghost && t->type != task_type_kick1 &&
         t->type != task_type_kick2 && t->type != task_type_timestep &&
-        t->type != task_type_drift)
+        t->type != task_type_drift_part && t->type != task_type_drift_gpart)
       continue;
 
     /* Get the task weight. */
@@ -557,7 +557,7 @@ static void repart_edge_metis(int partweights, int bothweights, int nodeID,
     /* Different weights for different tasks. */
     if (t->type == task_type_ghost || t->type == task_type_kick1 ||
         t->type == task_type_kick2 || t->type == task_type_timestep ||
-        t->type == task_type_drift) {
+        t->type == task_type_drift_part || t->type == task_type_drift_gpart) {
       /* Particle updates add only to vertex weight. */
       if (taskvweights) weights_v[cid] += w;
 

src/periodic.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016   Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PERIODIC_H
+#define SWIFT_PERIODIC_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Includes. */
+#include "inline.h"
+
+/**
+ * @brief Limits the value of x to be between a and b
+ *
+ * Only wraps once. If x > 2b, the returned value will be larger than b.
+ * Similarly for x < -b.
+ */
+#define box_wrap(x, a, b)                               \
+  ({                                                    \
+    const __typeof__(x) _x = (x);                       \
+    const __typeof__(a) _a = (a);                       \
+    const __typeof__(b) _b = (b);                       \
+    _x < _a ? (_x + _b) : ((_x > _b) ? (_x - _b) : _x); \
+  })
+
+/**
+ * @brief Find the smallest distance dx along one axis within a box of size
+ * box_size
+ *
+ * This macro evaluates its arguments exactly once.
+ *
+ * Only wraps once. If dx > 2b, the returned value will be larger than b.
+ * Similarly for dx < -b.
+ *
+ */
+__attribute__((always_inline)) INLINE static double nearest(double dx,
+                                                            double box_size) {
+  return dx > 0.5 * box_size ? (dx - box_size)
+                             : ((dx < -0.5 * box_size) ? (dx + box_size) : dx);
+}
+
+/**
+ * @brief Find the smallest distance dx along one axis within a box of size
+ * box_size
+ *
+ * This macro evaluates its arguments exactly once.
+ *
+ * Only wraps once. If dx > 2b, the returned value will be larger than b.
+ * Similarly for dx < -b.
+ *
+ */
+__attribute__((always_inline)) INLINE static float nearestf(float dx,
+                                                            float box_size) {
+  return dx > 0.5f * box_size
+             ? (dx - box_size)
+             : ((dx < -0.5f * box_size) ? (dx + box_size) : dx);
+}
+
+#endif /* SWIFT_PERIODIC_H */

src/runner.c

@@ -53,6 +53,7 @@
 #include "hydro_properties.h"
 #include "kick.h"
 #include "minmax.h"
+#include "runner_doiact_fft.h"
 #include "runner_doiact_vec.h"
 #include "scheduler.h"
 #include "sort_part.h"
@@ -331,7 +332,7 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
   TIMER_TIC;
 
   /* Check that the particles have been moved to the current time */
-  if (!cell_is_drifted(c, r->e)) error("Sorting un-drifted cell");
+  if (!cell_are_part_drifted(c, r->e)) error("Sorting un-drifted cell");
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Make sure the sort flags are consistent (downward). */
@@ -839,19 +840,35 @@ void runner_do_unskip_mapper(void *map_data, int num_elements,
   }
 }
 /**
- * @brief Drift particles in real space.
+ * @brief Drift all part in a cell.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_drift_part(struct runner *r, struct cell *c, int timer) {
+
+  TIMER_TIC;
+
+  cell_drift_part(c, r->e);
+
+  if (timer) TIMER_TOC(timer_drift_part);
+}
+
+/**
+ * @brief Drift all gpart in a cell.
  *
  * @param r The runner thread.
  * @param c The cell.
  * @param timer Are we timing this ?
  */
-void runner_do_drift_particles(struct runner *r, struct cell *c, int timer) {
+void runner_do_drift_gpart(struct runner *r, struct cell *c, int timer) {
 
   TIMER_TIC;
 
-  cell_drift_particles(c, r->e);
+  cell_drift_gpart(c, r->e);
 
-  if (timer) TIMER_TOC(timer_drift);
+  if (timer) TIMER_TOC(timer_drift_gpart);
 }
 
 /**
@@ -1522,7 +1539,7 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int clear_sorts,
   /* ... and store. */
   c->ti_end_min = ti_end_min;
   c->ti_end_max = ti_end_max;
-  c->ti_old = ti_current;
+  c->ti_old_part = ti_current;
   c->h_max = h_max;
 
   if (timer) TIMER_TOC(timer_dorecv_part);
@@ -1596,7 +1613,7 @@ void runner_do_recv_gpart(struct runner *r, struct cell *c, int timer) {
   /* ... and store. */
   c->ti_end_min = ti_end_min;
   c->ti_end_max = ti_end_max;
-  c->ti_old = ti_current;
+  c->ti_old_gpart = ti_current;
 
   if (timer) TIMER_TOC(timer_dorecv_gpart);
 
@@ -1669,7 +1686,7 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int timer) {
   /* ... and store. */
   c->ti_end_min = ti_end_min;
   c->ti_end_max = ti_end_max;
-  c->ti_old = ti_current;
+  c->ti_old_gpart = ti_current;
 
   if (timer) TIMER_TOC(timer_dorecv_spart);
 
@@ -1727,15 +1744,18 @@ 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 &&
             t->type != task_type_send && t->type != task_type_recv &&
-            t->type != task_type_kick1 && t->type != task_type_drift)
+            t->type != task_type_kick1 && t->type != task_type_drift_part &&
+            t->type != task_type_drift_gpart)
           error(
               "Task (type='%s/%s') should have been skipped ti_current=%lld "
               "c->ti_end_min=%lld",
@@ -1863,8 +1883,11 @@ void *runner_main(void *data) {
           runner_do_extra_ghost(r, ci, 1);
           break;
 #endif
-        case task_type_drift:
-          runner_do_drift_particles(r, ci, 1);
+        case task_type_drift_part:
+          runner_do_drift_part(r, ci, 1);
+          break;
+        case task_type_drift_gpart:
+          runner_do_drift_gpart(r, ci, 1);
           break;
         case task_type_kick1:
           runner_do_kick1(r, ci, 1);
@@ -1902,14 +1925,11 @@ void *runner_main(void *data) {
           }
           break;
 #endif
-        case task_type_grav_mm:
-          // runner_do_grav_mm(r, t->ci, 1);
-          break;
         case task_type_grav_down:
           runner_do_grav_down(r, t->ci, 1);
           break;
         case task_type_grav_top_level:
-          // runner_do_grav_top_level(r);
+          runner_do_grav_fft(r, 1);
           break;
         case task_type_grav_long_range:
           runner_do_grav_long_range(r, t->ci, 1);

src/runner.h

@@ -62,13 +62,15 @@ struct runner {
 void runner_do_ghost(struct runner *r, struct cell *c, int timer);
 void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer);
 void runner_do_sort(struct runner *r, struct cell *c, int flag, int clock);
-void runner_do_drift_particles(struct runner *r, struct cell *c, int timer);
+void runner_do_drift_part(struct runner *r, struct cell *c, int timer);
+void runner_do_drift_gpart(struct runner *r, struct cell *c, int timer);
 void runner_do_kick1(struct runner *r, struct cell *c, int timer);
 void runner_do_kick2(struct runner *r, struct cell *c, int timer);
 void runner_do_end_force(struct runner *r, struct cell *c, int timer);
 void runner_do_init(struct runner *r, struct cell *c, int timer);
 void runner_do_cooling(struct runner *r, struct cell *c, int timer);
 void runner_do_grav_external(struct runner *r, struct cell *c, int timer);
+void runner_do_grav_fft(struct runner *r, int timer);
 void *runner_main(void *data);
 void runner_do_unskip_mapper(void *map_data, int num_elements,
                              void *extra_data);

src/runner_doiact.h

@@ -903,7 +903,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
   /* Anything to do here? */
   if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
 
-  if (!cell_is_drifted(ci, e) || !cell_is_drifted(cj, e))
+  if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e))
     error("Interacting undrifted cells.");
 
   /* Get the sort ID. */
@@ -1147,7 +1147,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
   /* Anything to do here? */
   if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
 
-  if (!cell_is_drifted(ci, e) || !cell_is_drifted(cj, e))
+  if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e))
     error("Interacting undrifted cells.");
 
   /* Get the shift ID. */
@@ -1597,7 +1597,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
 
   if (!cell_is_active(c, e)) return;
 
-  if (!cell_is_drifted(c, e)) error("Interacting undrifted cell.");
+  if (!cell_are_part_drifted(c, e)) error("Interacting undrifted cell.");
 
   struct part *restrict parts = c->parts;
   const int count = c->count;
@@ -1846,7 +1846,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
 
   if (!cell_is_active(c, e)) return;
 
-  if (!cell_is_drifted(c, e)) error("Cell is not drifted");
+  if (!cell_are_part_drifted(c, e)) error("Cell is not drifted");
 
   struct part *restrict parts = c->parts;
   const int count = c->count;
@@ -2279,8 +2279,8 @@ void DOSUB_PAIR1(struct runner *r, struct cell *ci, struct cell *cj, int sid,
   else if (cell_is_active(ci, e) || cell_is_active(cj, e)) {
 
     /* Make sure both cells are drifted to the current timestep. */
-    if (!cell_is_drifted(ci, e)) cell_drift_particles(ci, e);
-    if (!cell_is_drifted(cj, e)) cell_drift_particles(cj, e);
+    if (!cell_are_part_drifted(ci, e)) cell_drift_part(ci, e);
+    if (!cell_are_part_drifted(cj, e)) cell_drift_part(cj, e);
 
     /* Do any of the cells need to be sorted first? */
     if (!(ci->sorted & (1 << sid)) ||
@@ -2333,7 +2333,7 @@ void DOSUB_SELF1(struct runner *r, struct cell *ci, int gettimer) {
   else {
 
     /* Drift the cell to the current timestep if needed. */
-    if (!cell_is_drifted(ci, r->e)) cell_drift_particles(ci, r->e);
+    if (!cell_are_part_drifted(ci, r->e)) cell_drift_part(ci, r->e);
 
 #if (DOSELF1 == runner_doself1_density) && defined(WITH_VECTORIZATION) && \
     defined(GADGET2_SPH)
@@ -2586,8 +2586,8 @@ void DOSUB_PAIR2(struct runner *r, struct cell *ci, struct cell *cj, int sid,
   else if (cell_is_active(ci, e) || cell_is_active(cj, e)) {
 
     /* Make sure both cells are drifted to the current timestep. */
-    if (!cell_is_drifted(ci, e)) cell_drift_particles(ci, e);
-    if (!cell_is_drifted(cj, e)) cell_drift_particles(cj, e);
+    if (!cell_are_part_drifted(ci, e)) cell_drift_part(ci, e);
+    if (!cell_are_part_drifted(cj, e)) cell_drift_part(cj, e);
 
     /* Do any of the cells need to be sorted first? */
     if (!(ci->sorted & (1 << sid)) ||
@@ -3218,7 +3218,7 @@ void DOSUB_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
       new_sid = sortlistID[new_sid];
 
       /* Do any of the cells need to be drifted first? */
-      if (!cell_is_drifted(cj, e)) cell_drift_particles(cj, e);
+      if (!cell_are_part_drifted(cj, e)) cell_drift_part(cj, e);
 
       DOPAIR_SUBSET(r, ci, parts, ind, count, cj);
     }

src/runner_doiact_fft.c

@@ -31,6 +31,283 @@
 #include "runner_doiact_fft.h"
 
 /* Local includes. */
+#include "engine.h"
+#include "error.h"
 #include "runner.h"
+#include "space.h"
+#include "timers.h"
 
-void runner_do_grav_fft(struct runner *r) {}
+/**
+ * @brief Returns 1D index of a 3D NxNxN array using row-major style.
+ *
+ * @param i Index along x.
+ * @param j Index along y.
+ * @param k Index along z.
+ * @param N Size of the array along one axis.
+ */
+__attribute__((always_inline)) INLINE static int row_major_id(int i, int j,
+                                                              int k, int N) {
+  return ((i % N) * N * N + (j % N) * N + (k % N));
+}
+
+/**
+ * @brief Assigns a given multipole to a density mesh using the CIC method.
+ *
+ * @param m The #multipole.
+ * @param rho The density mesh.
+ * @param N the size of the mesh along one axis.
+ * @param fac The width of a mesh cell.
+ */
+__attribute__((always_inline)) INLINE static void multipole_to_mesh_CIC(
+    const struct gravity_tensors* m, double* rho, int N, double fac) {
+
+  int i = (int)(fac * m->CoM[0]);
+  if (i >= N) i = N - 1;
+  const double dx = fac * m->CoM[0] - i;
+  const double tx = 1. - dx;
+
+  int j = (int)(fac * m->CoM[1]);
+  if (j >= N) j = N - 1;
+  const double dy = fac * m->CoM[1] - j;
+  const double ty = 1. - dy;
+
+  int k = (int)(fac * m->CoM[2]);
+  if (k >= N) k = N - 1;
+  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;
+  rho[row_major_id(i + 0, j + 1, k + 0, N)] += m->m_pole.M_000 * tx * dy * tz;
+  rho[row_major_id(i + 0, j + 1, k + 1, N)] += m->m_pole.M_000 * tx * dy * dz;
+  rho[row_major_id(i + 1, j + 0, k + 0, N)] += m->m_pole.M_000 * dx * ty * tz;
+  rho[row_major_id(i + 1, j + 0, k + 1, N)] += m->m_pole.M_000 * dx * ty * dz;
+  rho[row_major_id(i + 1, j + 1, k + 0, N)] += m->m_pole.M_000 * dx * dy * tz;
+  rho[row_major_id(i + 1, j + 1, k + 1, N)] += m->m_pole.M_000 * dx * dy * dz;
+}
+
+/**
+ * @brief Computes the potential on a multipole from a given mesh using the CIC
+ * method.
+ *
+ * @param m The #multipole.
+ * @param pot The potential mesh.
+ * @param N the size of the mesh along one axis.
+ * @param fac width of a mesh cell.
+ */
+__attribute__((always_inline)) INLINE static void mesh_to_multipole_CIC(
+    struct gravity_tensors* m, double* pot, int N, double fac) {
+
+  int i = (int)(fac * m->CoM[0]);
+  if (i >= N) i = N - 1;
+  const double dx = fac * m->CoM[0] - i;
+  const double tx = 1. - dx;
+
+  int j = (int)(fac * m->CoM[1]);
+  if (j >= N) j = N - 1;
+  const double dy = fac * m->CoM[1] - j;
+  const double ty = 1. - dy;
+
+  int k = (int)(fac * m->CoM[2]);
+  if (k >= N) k = N - 1;
+  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;
+  m->pot.F_000 += pot[row_major_id(i + 0, j + 1, k + 0, N)] * tx * dy * tz;
+  m->pot.F_000 += pot[row_major_id(i + 0, j + 1, k + 1, N)] * tx * dy * dz;
+  m->pot.F_000 += pot[row_major_id(i + 1, j + 0, k + 0, N)] * dx * ty * tz;
+  m->pot.F_000 += pot[row_major_id(i + 1, j + 0, k + 1, N)] * dx * ty * dz;
+  m->pot.F_000 += pot[row_major_id(i + 1, j + 1, k + 0, N)] * dx * dy * tz;
+  m->pot.F_000 += pot[row_major_id(i + 1, j + 1, k + 1, N)] * dx * dy * dz;
+}
+
+/**
+ * @brief Computes the potential on the top multipoles using a Fourier transform
+ *
+ * @param r The #runner task
+ * @param timer Are we timing this ?
+ */
+void runner_do_grav_fft(struct runner* r, int timer) {
+
+#ifdef HAVE_FFTW
+
+  const struct engine* e = r->e;
+  const struct space* s = e->s;
+  const integertime_t ti_current = e->ti_current;
+  const double a_smooth = e->gravity_properties->a_smooth;
+  const double box_size = s->dim[0];
+  const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
+
+  TIMER_TIC;
+
+  if (cdim[0] != cdim[1] || cdim[0] != cdim[2]) error("Non-square mesh");
+
+  /* Some useful constants */
+  const int N = cdim[0];
+  const int N_half = N / 2;
+  const double cell_fac = N / box_size;
+
+  /* Recover the list of top-level multipoles */