Fix to #363. Drifting decision needs to use the same multipole information as...

Fix to #363. Drifting decision needs to use the same multipole information as the actual tree walk. Also fix some gravity_cache padding values.

examples/EAGLE_6/eagle_6.yml

@@ -13,9 +13,6 @@ TimeIntegration:
   dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-4  # The maximal time-step size of the simulation (in internal units).
 
-Scheduler:
-  cell_split_size: 64
-  
 # Parameters governing the snapshots
 Snapshots:
   basename:            eagle # Common part of the name of output files
@@ -29,8 +26,8 @@ Statistics:
 # Parameters for the self-gravity scheme
 Gravity:
   eta:                   0.025    # Constant dimensionless multiplier for time integration.
-  theta:                 0.7      # Opening angle (Multipole acceptance criterion)
-  epsilon:               0.001   # Softening length (in internal units).
+  theta:                 0.85     # Opening angle (Multipole acceptance criterion)
+  epsilon:               0.001    # Softening length (in internal units).
   
 # Parameters for the hydrodynamics scheme
 SPH:

src/cell.c

@@ -1692,6 +1692,19 @@ 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) {
 

src/engine.c

@@ -3656,7 +3656,7 @@ void engine_step(struct engine *e) {
   if (e->nodeID == 0) {
 
     /* Print some information to the screen */
-    printf("  %6d %14e %14e %10zu %10zu %10zu %21.3f\n", e->step, e->time,
+    printf("  %6d %10lld %3d %14e %14e %10zu %10zu %10zu %21.3f\n", e->step, e->ti_current, e->max_active_bin, e->time,
            e->timeStep, e->updates, e->g_updates, e->s_updates,
            e->wallclock_time);
     fflush(stdout);

src/gravity_cache.h

@@ -189,13 +189,14 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
    * that have a reasonable magnitude. We use the cell width for this */
   const float pos_padded[3] = {-2. * cell->width[0], -2. * cell->width[1],
                                -2. * cell->width[2]};
+  const float eps_padded = epsilon[0];
 
   /* Pad the caches */
   for (int i = gcount; i < gcount_padded; ++i) {
     x[i] = pos_padded[0];
     y[i] = pos_padded[1];
     z[i] = pos_padded[2];
-    epsilon[i] = 0.f;
+    epsilon[i] = eps_padded;
     m[i] = 0.f;
     active[i] = 0;
     use_mpole[i] = 0;
@@ -248,12 +249,14 @@ gravity_cache_populate_no_mpole(timebin_t max_active_bin,
    * that have a reasonable magnitude. We use the cell width for this */
   const float pos_padded[3] = {-2. * cell->width[0], -2. * cell->width[1],
                                -2. * cell->width[2]};
+  const float eps_padded = epsilon[0];
+
   /* Pad the caches */
   for (int i = gcount; i < gcount_padded; ++i) {
     x[i] = pos_padded[0];
     y[i] = pos_padded[1];
     z[i] = pos_padded[2];
-    epsilon[i] = 0.f;
+    epsilon[i] = eps_padded;
     m[i] = 0.f;
     active[i] = 0;
   }

src/multipole.h

@@ -188,12 +188,18 @@ 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;
     };
@@ -232,7 +238,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 + 2. * x_diff;
+  m->r_max = m->r_max_rebuild + 0. * x_diff;
 }
 
 /**

src/runner_doiact_grav.h

@@ -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[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];
+  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];
 
   /* Apply BC */
   if (periodic) {
@@ -999,7 +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, multi_j->r_max, theta_crit2, r2)) {
+  if (gravity_M2L_accept(multi_i->r_max_old, multi_j->r_max_old, theta_crit2, r2)) {
 
     /* MATTHIEU: make a symmetric M-M interaction function ! */
     runner_dopair_grav_mm(r, ci, cj);
@@ -1012,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;
-    const double rj_max = multi_j->r_max;
+    const double ri_max = multi_i->r_max_old;
+    const double rj_max = multi_j->r_max_old;
 
     /* Split the larger of the two cells and start over again */
     if (ri_max > rj_max) {
@@ -1169,7 +1169,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[0], multi_i->CoM[1], multi_i->CoM[2]};
+  const double CoM_i[3] = {multi_i->CoM_old[0], multi_i->CoM_old[1], multi_i->CoM_old[2]};
   const double CoM_rebuild_i[3] = {multi_i->CoM_rebuild[0],
                                    multi_i->CoM_rebuild[1],
                                    multi_i->CoM_rebuild[2]};
@@ -1186,9 +1186,9 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
     if (ci == cj || cj->gcount == 0) continue;
 
     /* Get the distance between the CoMs */
-    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];
+    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];
 
     /* Apply BC */
     if (periodic) {
@@ -1209,7 +1209,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
     }
 
     /* Check the multipole acceptance criterion */
-    if (gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2)) {
+    if (gravity_M2L_accept(multi_i->r_max_old, multi_j->r_max_old, theta_crit2, r2)) {
 
       /* Go for a (non-symmetric) M-M calculation */
       runner_dopair_grav_mm(r, ci, cj);

src/vector.h

@@ -433,7 +433,7 @@ __attribute__((always_inline)) INLINE vector vec_reciprocal_sqrt(vector x) {
 
 #else
 /* Needed for cache alignment. */
-#define VEC_SIZE 16
+#define VEC_SIZE 8
 #endif /* WITH_VECTORIZATION */
 
 #endif /* VEC_MACRO */