Drift multipoles before using them for interactions.

src/cell.c

@@ -1503,7 +1503,21 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
  * @param e The #engine (to get ti_current).
  */
 void cell_drift_multipole(struct cell *c, const struct engine *e) {
-  error("To be implemented");
+
+  const double timeBase = e->timeBase;
+  const integertime_t ti_old_multipole = c->ti_old_multipole;
+  const integertime_t ti_current = e->ti_current;
+
+  /* Drift from the last time the cell was drifted to the current time */
+  const double dt = (ti_current - ti_old_multipole) * timeBase;
+
+  /* Check that we are actually going to move forward. */
+  if (ti_current < ti_old_multipole) error("Attempt to drift to the past");
+
+  if (ti_current > ti_old_multipole) gravity_drift(c->multipole, dt);
+
+  /* Update the time of the last drift */
+  c->ti_old_multipole = ti_current;
 }
 
 /**

src/gravity/Default/gravity_iact.h

@@ -44,6 +44,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
   const float u = r * rlr_inv;
   float f_lr, fi, fj, W;
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (r == 0.f) error("Interacting particles with 0 distance");
+#endif
+
   /* Get long-range correction */
   kernel_long_grav_eval(u, &f_lr);
 
@@ -107,6 +111,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
   const float u = r * rlr_inv;
   float f_lr, f, W;
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (r == 0.f) error("Interacting particles with 0 distance");
+#endif
+
   /* Get long-range correction */
   kernel_long_grav_eval(u, &f_lr);
 
@@ -142,52 +150,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
     const struct multipole *multi) {
 
   error("Dead function");
-
-  /*   /\* Apply the gravitational acceleration. *\/ */
-  /*   const float r = sqrtf(r2); */
-  /*   const float ir = 1.f / r; */
-  /*   const float u = r * rlr_inv; */
-  /*   const float mrinv3 = multi->mass * ir * ir * ir; */
-
-  /*   /\* Get long-range correction *\/ */
-  /*   float f_lr; */
-  /*   kernel_long_grav_eval(u, &f_lr); */
-
-  /* #if const_gravity_multipole_order < 2 */
-
-  /*   /\* 0th and 1st order terms *\/ */
-  /*   gp->a_grav[0] += mrinv3 * f_lr * dx[0]; */
-  /*   gp->a_grav[1] += mrinv3 * f_lr * dx[1]; */
-  /*   gp->a_grav[2] += mrinv3 * f_lr * dx[2]; */
-
-  /* #elif const_gravity_multipole_order == 2 */
-  /*   /\* Terms up to 2nd order (quadrupole) *\/ */
-
-  /*   /\* Follows the notation in Bonsai *\/ */
-  /*   const float mrinv5 = mrinv3 * ir * ir; */
-  /*   const float mrinv7 = mrinv5 * ir * ir; */
-
-  /*   const float D1 = -mrinv3; */
-  /*   const float D2 = 3.f * mrinv5; */
-  /*   const float D3 = -15.f * mrinv7; */
-
-  /*   const float q = multi->I_xx + multi->I_yy + multi->I_zz; */
-  /*   const float qRx = */
-  /*       multi->I_xx * dx[0] + multi->I_xy * dx[1] + multi->I_xz * dx[2]; */
-  /*   const float qRy = */
-  /*       multi->I_xy * dx[0] + multi->I_yy * dx[1] + multi->I_yz * dx[2]; */
-  /*   const float qRz = */
-  /*       multi->I_xz * dx[0] + multi->I_yz * dx[1] + multi->I_zz * dx[2]; */
-  /*   const float qRR = qRx * dx[0] + qRy * dx[1] + qRz * dx[2]; */
-  /*   const float C = D1 + 0.5f * D2 * q + 0.5f * D3 * qRR; */
-
-  /*   gp->a_grav[0] -= f_lr * (C * dx[0] + D2 * qRx); */
-  /*   gp->a_grav[1] -= f_lr * (C * dx[1] + D2 * qRy); */
-  /*   gp->a_grav[2] -= f_lr * (C * dx[2] + D2 * qRz); */
-
-  /* #else */
-  /* #error "Multipoles of order >2 not yet implemented." */
-  /* #endif */
 }
 
 #endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */

src/multipole.h

@@ -1127,6 +1127,11 @@ INLINE static void gravity_L2L(struct grav_tensor *la,
   /* Initialise everything to zero */
   gravity_field_tensors_init(la);
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (lb->num_interacted == 0) error("Shifting tensors that did not interact");
+  la->num_interacted = lb->num_interacted;
+#endif
+
   /* Distance to shift by */
   const double dx[3] = {pos_a[0] - pos_b[0], pos_a[1] - pos_b[1],
                         pos_a[2] - pos_b[2]};
@@ -1222,11 +1227,6 @@ INLINE static void gravity_L2L(struct grav_tensor *la,
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 3
 #error "Missing implementation for order >3"
 #endif
-
-#ifdef SWIFT_DEBUG_CHECKS
-  if (lb->num_interacted == 0) error("Shifting tensors that did not interact");
-  la->num_interacted = lb->num_interacted;
-#endif
 }
 
 /**

src/runner.c

@@ -1375,8 +1375,8 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
 #ifdef SWIFT_DEBUG_CHECKS
       if (e->policy & engine_policy_self_gravity) {
         gp->num_interacted++;
-        if (gp->num_interacted - e->s->nr_gparts)
-          error(
+        if (gp->num_interacted != (long long)e->s->nr_gparts)
+          message(
               "g-particle did not interact gravitationally with all other "
               "particles gp->num_interacted=%lld, total_gparts=%zd",
               gp->num_interacted, e->s->nr_gparts);

src/runner_doiact_grav.h

@@ -84,9 +84,8 @@ void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
  * @param ci The #cell with field tensor to interact.
  * @param cj The #cell with the multipole.
  */
-void runner_dopair_grav_mm(const struct runner *r,
-                           const struct cell *restrict ci,
-                           const struct cell *restrict cj) {
+void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
+                           struct cell *restrict cj) {
 
   const struct engine *e = r->e;
   const int periodic = e->s->periodic;
@@ -97,12 +96,18 @@ void runner_dopair_grav_mm(const struct runner *r,
   TIMER_TIC;
 
 #ifdef SWIFT_DEBUG_CHECKS
+  if (ci == cj) error("Interacting a cell with itself using M2L");
+
   if (multi_j->M_000 == 0.f) error("Multipole does not seem to have been set.");
 #endif
 
   /* Anything to do here? */
   if (!cell_is_active(ci, e)) return;
 
+  /* Do we need to drift the multipole ? */
+  if (cj->ti_old_multipole != e->ti_current) cell_drift_multipole(cj, e);
+
+  /* Let's interact at this level */
   gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
               cj->multipole->CoM, periodic * 0);
 
@@ -147,7 +152,7 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
   TIMER_TIC;
 
 #ifdef SWIFT_DEBUG_CHECKS
-  if (ci->width[0] != cj->width[0])  // MATTHIEU sanity check
+  if (ci->width[0] != cj->width[0])
     error("Non matching cell sizes !! h_i=%f h_j=%f", ci->width[0],
           cj->width[0]);
 #endif
@@ -533,6 +538,9 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
   /* Anything to do here? */
   if (!cell_is_active(ci, e)) return;
 
+  /* Drift our own multipole if need be */
+  if (ci->ti_old_multipole != e->ti_current) cell_drift_multipole(ci, e);
+
   /* Loop over all the cells and go for a p-m interaction if far enough but not
    * too far */
   for (int i = 0; i < nr_cells; ++i) {