Merge branch 'gravity_fixes' into 'master'

Fixes to the gravity code

This fixes #363. The long-range gravity task has to consider all the cells that were assigned to it at the rebuild time. 

Additional improvements:
 - Drift the multipoles while activating the tasks,
 - Avoid FPEs in the FFT task,
 - Adjust the number of tasks per top-level required to run.

See merge request !430

src/cell.c

@@ -1692,18 +1692,6 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
   const double dim[3] = {sp->dim[0], sp->dim[1], sp->dim[2]};
   const double theta_crit2 = e->gravity_properties->theta_crit2;
 
-  /* Store the current dx_max and h_max values. */
-  ci->multipole->r_max_old = ci->multipole->r_max;
-  ci->multipole->CoM_old[0] = ci->multipole->CoM[0];
-  ci->multipole->CoM_old[1] = ci->multipole->CoM[1];
-  ci->multipole->CoM_old[2] = ci->multipole->CoM[2];
-  if (cj != NULL) {
-    cj->multipole->r_max_old = cj->multipole->r_max;
-    cj->multipole->CoM_old[0] = cj->multipole->CoM[0];
-    cj->multipole->CoM_old[1] = cj->multipole->CoM[1];
-    cj->multipole->CoM_old[2] = cj->multipole->CoM[2];
-  }
-
   /* Self interaction? */
   if (cj == NULL) {
 
@@ -1736,6 +1724,9 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
     /* Anything to do here? */
     if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
 
+    if (ci->ti_old_multipole < e->ti_current) cell_drift_multipole(ci, e);
+    if (cj->ti_old_multipole < e->ti_current) cell_drift_multipole(cj, e);
+
     /* Recover the multipole information */
     struct gravity_tensors *const multi_i = ci->multipole;
     struct gravity_tensors *const multi_j = cj->multipole;

src/engine.c

@@ -2971,7 +2971,7 @@ int engine_estimate_nr_tasks(struct engine *e) {
 #endif
   }
   if (e->policy & engine_policy_self_gravity) {
-    n1 += 24;
+    n1 += 32;
     n2 += 1;
 #ifdef WITH_MPI
     n2 += 2;
@@ -3070,13 +3070,6 @@ void engine_rebuild(struct engine *e, int clean_h_values) {
   /* Re-build the tasks. */
   engine_maketasks(e);
 
-  /* Run through the tasks and mark as skip or not. */
-  if (engine_marktasks(e))
-    error("engine_marktasks failed after space_rebuild.");
-
-/* Print the status of the system */
-// if (e->verbose) engine_print_task_counts(e);
-
 #ifdef SWIFT_DEBUG_CHECKS
   /* Check that all cells have been drifted to the current time.
    * That can include cells that have not
@@ -3085,6 +3078,13 @@ void engine_rebuild(struct engine *e, int clean_h_values) {
                           e->policy & engine_policy_self_gravity);
 #endif
 
+  /* Run through the tasks and mark as skip or not. */
+  if (engine_marktasks(e))
+    error("engine_marktasks failed after space_rebuild.");
+
+  /* Print the status of the system */
+  // if (e->verbose) engine_print_task_counts(e);
+
   if (e->verbose)
     message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
             clocks_getunit());
@@ -3692,8 +3692,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 ? */
@@ -3901,7 +3901,7 @@ void engine_do_drift_top_multipoles_mapper(void *map_data, int num_elements,
     if (c != NULL && c->nodeID == e->nodeID) {
 
       /* Drift the multipole at this level only */
-      cell_drift_multipole(c, e);
+      if (c->ti_old_multipole != e->ti_current) cell_drift_multipole(c, e);
     }
   }
 }

src/gravity/Default/gravity_debug.h

@@ -26,6 +26,10 @@ __attribute__((always_inline)) INLINE static void gravity_debug_particle(
       "x=[%.5e,%.5e,%.5e], v_full=[%.5e,%.5e,%.5e], a=[%.5e,%.5e,%.5e]\n",
       p->mass, p->epsilon, p->time_bin, p->x[0], p->x[1], p->x[2], p->v_full[0],
       p->v_full[1], p->v_full[2], p->a_grav[0], p->a_grav[1], p->a_grav[2]);
+#ifdef SWIFT_DEBUG_CHECKS
+  printf("num_interacted=%lld ti_drift=%lld ti_kick=%lld\n", p->num_interacted,
+         p->ti_drift, p->ti_kick);
+#endif
 }
 
 #endif /* SWIFT_DEFAULT_GRAVITY_DEBUG_H */

src/kernel_long_gravity.h

@@ -78,7 +78,7 @@ __attribute__((always_inline)) INLINE static void fourier_kernel_long_grav_eval(
 #else
   const double u = sqrt(u2);
   const double arg = M_PI_2 * u;
-  *W = arg / sinh(arg);
+  *W = arg / (sinh(arg) + FLT_MIN);
 #endif
 }
 

src/multipole.h

@@ -188,18 +188,12 @@ struct gravity_tensors {
       /*! Centre of mass of the matter dsitribution */
       double CoM[3];
 
-      /*! Centre of mass of the matter dsitribution at the last time-step */
-      double CoM_old[3];
-
       /*! Centre of mass of the matter dsitribution at the last rebuild */
       double CoM_rebuild[3];
 
       /*! Upper limit of the CoM<->gpart distance */
       double r_max;
 
-      /*! Upper limit of the CoM<->gpart distance at the last time-step */
-      double r_max_old;
-
       /*! Upper limit of the CoM<->gpart distance at the last rebuild */
       double r_max_rebuild;
     };
@@ -238,7 +232,7 @@ INLINE static void gravity_drift(struct gravity_tensors *m, double dt,
   m->CoM[2] += dz;
 
   /* Conservative change in maximal radius containing all gpart */
-  m->r_max = m->r_max_rebuild + 0. * x_diff;
+  m->r_max = m->r_max_rebuild + x_diff;
 }
 
 /**

src/runner_doiact_fft.c

@@ -231,7 +231,7 @@ void runner_do_grav_fft(struct runner* r, int timer) {
         const int kz = (k > N_half ? k - N : k);
         const double kz_d = (double)kz;
         const double fz = k_fac * kz_d;
-        const double sinc_kz_inv = (kz != 0) ? fz / sin(fz) : 1.;
+        const double sinc_kz_inv = (kz != 0) ? fz / (sin(fz) + FLT_MIN) : 1.;
 
         /* Norm of vector in Fourier space */
         const double k2 = (kx_d * kx_d + ky_d * ky_d + kz_d * kz_d);
@@ -242,7 +242,7 @@ void runner_do_grav_fft(struct runner* r, int timer) {
         /* Green function */
         double W;
         fourier_kernel_long_grav_eval(k2 * a_smooth2, &W);
-        const double green_cor = green_fac * W / k2;
+        const double green_cor = green_fac * W / (k2 + FLT_MIN);
 
         /* Deconvolution of CIC */
         const double CIC_cor = sinc_kx_inv * sinc_ky_inv * sinc_kz_inv;

src/runner_doiact_grav.h

@@ -155,7 +155,7 @@ void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
 #endif
 
   /* Do we need to drift the multipole ? */
-  if (cj->ti_old_multipole != e->ti_current) cell_drift_multipole(cj, e);
+  if (cj->ti_old_multipole != e->ti_current) error("Undrifted multipole");
 
   /* Let's interact at this level */
   gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
@@ -469,9 +469,9 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
 
   /* Do we need to drift the multipoles ? */
   if (cj_active && ci->ti_old_multipole != e->ti_current)
-    cell_drift_multipole(ci, e);
+    error("Un-drifted multipole");
   if (ci_active && cj->ti_old_multipole != e->ti_current)
-    cell_drift_multipole(cj, e);
+    error("Un-drifted multipole");
 
   /* Centre of the cell pair */
   const double loc[3] = {ci->loc[0],   // + 0. * ci->width[0],
@@ -970,9 +970,9 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
   struct gravity_tensors *const multi_j = cj->multipole;
 
   /* Get the distance between the CoMs */
-  double dx = multi_i->CoM_old[0] - multi_j->CoM_old[0];
-  double dy = multi_i->CoM_old[1] - multi_j->CoM_old[1];
-  double dz = multi_i->CoM_old[2] - multi_j->CoM_old[2];
+  double dx = multi_i->CoM[0] - multi_j->CoM[0];
+  double dy = multi_i->CoM[1] - multi_j->CoM[1];
+  double dz = multi_i->CoM[2] - multi_j->CoM[2];
 
   /* Apply BC */
   if (periodic) {
@@ -999,8 +999,7 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
    * option... */
 
   /* Can we use M-M interactions ? */
-  if (gravity_M2L_accept(multi_i->r_max_old, multi_j->r_max_old, theta_crit2,
-                         r2)) {
+  if (gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2)) {
 
     /* MATTHIEU: make a symmetric M-M interaction function ! */
     runner_dopair_grav_mm(r, ci, cj);
@@ -1013,8 +1012,8 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
   /* Alright, we'll have to split and recurse. */
   else {
 
-    const double ri_max = multi_i->r_max_old;
-    const double rj_max = multi_j->r_max_old;
+    const double ri_max = multi_i->r_max;
+    const double rj_max = multi_j->r_max;
 
     /* Split the larger of the two cells and start over again */
     if (ri_max > rj_max) {
@@ -1131,19 +1130,6 @@ void runner_doself_grav(struct runner *r, struct cell *c, int gettimer) {
  */
 void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
 
-#if ICHECK > 0
-  for (int pid = 0; pid < ci->gcount; pid++) {
-
-    /* Get a hold of the ith part in ci. */
-    struct gpart *restrict gp = &ci->gparts[pid];
-
-    if (gp->id_or_neg_offset == ICHECK)
-      message("id=%lld loc=[ %f %f %f ] size= %f count= %d",
-              gp->id_or_neg_offset, ci->loc[0], ci->loc[1], ci->loc[2],
-              ci->width[0], ci->gcount);
-  }
-#endif
-
   /* Some constants */
   const struct engine *e = r->e;
   const struct space *s = e->s;
@@ -1170,8 +1156,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
 
   /* Recover the local multipole */
   struct gravity_tensors *const multi_i = ci->multipole;
-  const double CoM_i[3] = {multi_i->CoM_old[0], multi_i->CoM_old[1],
-                           multi_i->CoM_old[2]};
+  const double CoM_i[3] = {multi_i->CoM[0], multi_i->CoM[1], multi_i->CoM[2]};
   const double CoM_rebuild_i[3] = {multi_i->CoM_rebuild[0],
                                    multi_i->CoM_rebuild[1],
                                    multi_i->CoM_rebuild[2]};
@@ -1187,65 +1172,58 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
     /* Avoid stupid cases */
     if (ci == cj || cj->gcount == 0) continue;
 
-    /* Get the distance between the CoMs */
-    double dx = CoM_i[0] - multi_j->CoM_old[0];
-    double dy = CoM_i[1] - multi_j->CoM_old[1];
-    double dz = CoM_i[2] - multi_j->CoM_old[2];
+    /* Get the distance between the CoMs at the last rebuild*/
+    double dx_r = CoM_rebuild_i[0] - multi_j->CoM_rebuild[0];
+    double dy_r = CoM_rebuild_i[1] - multi_j->CoM_rebuild[1];
+    double dz_r = CoM_rebuild_i[2] - multi_j->CoM_rebuild[2];
 
     /* Apply BC */
     if (periodic) {
-      dx = nearest(dx, dim[0]);
-      dy = nearest(dy, dim[1]);
-      dz = nearest(dz, dim[2]);
-    }
-    const double r2 = dx * dx + dy * dy + dz * dz;
-
-    /* Are we beyond the distance where the truncated forces are 0 ?*/
-    if (periodic && r2 > max_distance2) {
-
-#ifdef SWIFT_DEBUG_CHECKS
-      /* Need to account for the interactions we missed */
-      multi_i->pot.num_interacted += multi_j->m_pole.num_gpart;
-#endif
-      continue;
+      dx_r = nearest(dx_r, dim[0]);
+      dy_r = nearest(dy_r, dim[1]);
+      dz_r = nearest(dz_r, dim[2]);
     }
+    const double r2_rebuild = dx_r * dx_r + dy_r * dy_r + dz_r * dz_r;
 
-    /* Check the multipole acceptance criterion */
-    if (gravity_M2L_accept(multi_i->r_max_old, multi_j->r_max_old, theta_crit2,
-                           r2)) {
+    /* Are we in charge of this cell pair? */
+    if (gravity_M2L_accept(multi_i->r_max_rebuild, multi_j->r_max_rebuild,
+                           theta_crit2, r2_rebuild)) {
 
-      /* Go for a (non-symmetric) M-M calculation */
-      runner_dopair_grav_mm(r, ci, cj);
-
-    } else {
-
-      /* Let's check whether we need to still operate on this pair */
-
-      /* Get the distance between the CoMs at the last rebuild*/
-      double dx_rebuild = CoM_rebuild_i[0] - multi_j->CoM_rebuild[0];
-      double dy_rebuild = CoM_rebuild_i[1] - multi_j->CoM_rebuild[1];
-      double dz_rebuild = CoM_rebuild_i[2] - multi_j->CoM_rebuild[2];
+      /* Let's compute the current distance between the cell pair*/
+      double dx = CoM_i[0] - multi_j->CoM[0];
+      double dy = CoM_i[1] - multi_j->CoM[1];
+      double dz = CoM_i[2] - multi_j->CoM[2];
 
       /* Apply BC */
       if (periodic) {
-        dx_rebuild = nearest(dx_rebuild, dim[0]);
-        dy_rebuild = nearest(dy_rebuild, dim[1]);
-        dz_rebuild = nearest(dz_rebuild, dim[2]);
+        dx = nearest(dx, dim[0]);
+        dy = nearest(dy, dim[1]);
+        dz = nearest(dz, dim[2]);
       }
-      const double r2_rebuild = dx_rebuild * dx_rebuild +
-                                dy_rebuild * dy_rebuild +
-                                dz_rebuild * dz_rebuild;
+      const double r2 = dx * dx + dy * dy + dz * dz;
 
-      /* Is the criterion violated now but was OK at the last rebuild ? */
-      if (gravity_M2L_accept(multi_i->r_max_rebuild, multi_j->r_max_rebuild,
-                             theta_crit2, r2_rebuild)) {
+      /* Are we beyond the distance where the truncated forces are 0 ?*/
+      if (periodic && r2 > max_distance2) {
 
+#ifdef SWIFT_DEBUG_CHECKS
+        /* Need to account for the interactions we missed */
+        multi_i->pot.num_interacted += multi_j->m_pole.num_gpart;
+#endif
+        continue;
+      }
+
+      /* Check the multipole acceptance criterion */
+      if (gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2)) {
+
+        /* Go for a (non-symmetric) M-M calculation */
+        runner_dopair_grav_mm(r, ci, cj);
+      } else {
         /* Alright, we have to take charge of that pair in a different way. */
         // MATTHIEU: We should actually open the tree-node here and recurse.
         runner_dopair_grav_mm(r, ci, cj);
       }
-    }
-  } /* Loop over top-level cells */
+    } /* We are in charge of this pair */
+  }   /* Loop over top-level cells */
 
   if (timer) TIMER_TOC(timer_dograv_long_range);
 }