Merge branch 'fix_pair_vec' into 'master'

 - Use the correct value of 'h_max' in the construction of the caches in the vectorized DOPAIR1 routine. (fix #312)
 - Update the test27 unit test to also randomly perturb the smoothing length values.
 - Set a correct value for h_max in the test27 case and not just a loose upper-bound.
 - Call the non-vectorized version of DOPAIR1 for the corner interactions.

See merge request !348

EAGLE_25 now produces the same number of time-steps whether running with or without vectorization. 

src/runner.c

@@ -1818,13 +1818,8 @@ void *runner_main(void *data) {
           break;
 
         case task_type_pair:
-          if (t->subtype == task_subtype_density) {
-#if defined(WITH_VECTORIZATION) && defined(GADGET2_SPH)
-            runner_dopair1_density_vec(r, ci, cj);
-#else
-            runner_dopair1_density(r, ci, cj);
-#endif
-          }
+          if (t->subtype == task_subtype_density)
+            runner_dopair1_branch_density(r, ci, cj);
 #ifdef EXTRA_HYDRO_LOOP
           else if (t->subtype == task_subtype_gradient)
             runner_dopair1_gradient(r, ci, cj);

src/runner_doiact.h

@@ -884,8 +884,11 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci,
  * @param r The #runner.
  * @param ci The first #cell.
  * @param cj The second #cell.
+ * @param sid The direction of the pair
+ * @param shift The shift vector to apply to the particles in ci.
  */
-void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
+void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
+             const double *shift) {
 
   const struct engine *restrict e = r->e;
 
@@ -900,22 +903,6 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
 
   TIMER_TIC;
 
-  /* Anything to do here? */
-  if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
-
-  if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e))
-    error("Interacting undrifted cells.");
-
-  /* Get the sort ID. */
-  double shift[3] = {0.0, 0.0, 0.0};
-  const int sid = space_getsid(e->s, &ci, &cj, shift);
-
-  /* Have the cells been sorted? */
-  if (!(ci->sorted & (1 << sid)) || ci->dx_max_sort > space_maxreldx * ci->dmin)
-    runner_do_sort(r, ci, (1 << sid), 1);
-  if (!(cj->sorted & (1 << sid)) || cj->dx_max_sort > space_maxreldx * cj->dmin)
-    runner_do_sort(r, cj, (1 << sid), 1);
-
   /* Get the cutoff shift. */
   double rshift = 0.0;
   for (int k = 0; k < 3; k++) rshift += shift[k] * runner_shift[sid][k];
@@ -1116,6 +1103,51 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
   TIMER_TOC(TIMER_DOPAIR);
 }
 
+/**
+ * @brief Determine which version of DOPAIR1 needs to be called depending on the
+ * orientation of the cells or whether DOPAIR1 needs to be called at all.
+ *
+ * @param r #runner
+ * @param ci #cell ci
+ * @param cj #cell cj
+ *
+ */
+void DOPAIR1_BRANCH(struct runner *r, struct cell *ci, struct cell *cj) {
+
+  const struct engine *restrict e = r->e;
+
+  /* Anything to do here? */
+  if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
+
+  /* Check that cells are drifted. */
+  if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e))
+    error("Interacting undrifted cells.");
+
+  /* Get the sort ID. */
+  double shift[3] = {0.0, 0.0, 0.0};
+  const int sid = space_getsid(e->s, &ci, &cj, shift);
+
+  /* Have the cells been sorted? */
+  if (!(ci->sorted & (1 << sid)) || ci->dx_max_sort > space_maxreldx * ci->dmin)
+    runner_do_sort(r, ci, (1 << sid), 1);
+  if (!(cj->sorted & (1 << sid)) || cj->dx_max_sort > space_maxreldx * cj->dmin)
+    runner_do_sort(r, cj, (1 << sid), 1);
+
+  /* Have the cells been sorted? */
+  if (!(ci->sorted & (1 << sid)) || !(cj->sorted & (1 << sid)))
+    error("Trying to interact unsorted cells.");
+
+#if defined(WITH_VECTORIZATION) && defined(GADGET2_SPH) && \
+    (DOPAIR1_BRANCH == runner_dopair1_density_branch)
+  if (!sort_is_corner(sid))
+    runner_dopair1_density_vec(r, ci, cj, sid, shift);
+  else
+    DOPAIR1(r, ci, cj, sid, shift);
+#else
+  DOPAIR1(r, ci, cj, sid, shift);
+#endif
+}
+
 /**
  * @brief Compute the interactions between a cell pair (symmetric)
  *
@@ -2290,13 +2322,8 @@ void DOSUB_PAIR1(struct runner *r, struct cell *ci, struct cell *cj, int sid,
         cj->dx_max_sort > cj->dmin * space_maxreldx)
       runner_do_sort(r, cj, (1 << sid), 1);
 
-/* Compute the interactions. */
-#if (DOPAIR1 == runner_dopair1_density) && defined(WITH_VECTORIZATION) && \
-    defined(GADGET2_SPH)
-    runner_dopair1_density_vec(r, ci, cj);
-#else
-    DOPAIR1(r, ci, cj);
-#endif
+    /* Compute the interactions. */
+    DOPAIR1_BRANCH(r, ci, cj);
   }
 
   if (gettimer) TIMER_TOC(TIMER_DOSUB_PAIR);

src/runner_doiact_vec.c

@@ -20,13 +20,12 @@
 /* Config parameters. */
 #include "../config.h"
 
-#include "swift.h"
-
-#include "active.h"
-
 /* This object's header. */
 #include "runner_doiact_vec.h"
 
+/* Local headers. */
+#include "active.h"
+
 #ifdef WITH_VECTORIZATION
 /**
  * @brief Compute the vector remainder interactions from the secondary cache.
@@ -262,40 +261,40 @@ __attribute__((always_inline)) INLINE static void storeInteractions(
   }
 }
 
-/* @brief Populates the arrays max_di and max_dj with the maximum distances of
+/**
+ * @brief Populates the arrays max_di and max_dj with the maximum distances of
  * particles into their neighbouring cells. Also finds the first pi that
  * interacts with any particle in cj and the last pj that interacts with any
  * particle in ci.
+ *
  * @param ci #cell pointer to ci
  * @param cj #cell pointer to cj
  * @param sort_i #entry array for particle distance in ci
  * @param sort_j #entry array for particle distance in cj
- * @param ci_cache #cache for cell ci
- * @param cj_cache #cache for cell cj
  * @param dx_max maximum particle movement allowed in cell
  * @param rshift cutoff shift
+ * @param hi_max Maximal smoothing length in cell ci
+ * @param hj_max Maximal smoothing length in cell cj
+ * @param di_max Maximal position on the axis that can interact in cell ci
+ * @param dj_min Minimal position on the axis that can interact in cell ci
  * @param max_di array to hold the maximum distances of pi particles into cell
  * cj
  * @param max_dj array to hold the maximum distances of pj particles into cell
  * cj
  * @param init_pi first pi to interact with a pj particle
  * @param init_pj last pj to interact with a pi particle
+ * @param e The #engine.
  */
 __attribute__((always_inline)) INLINE static void populate_max_d_no_cache(
     const struct cell *ci, const struct cell *cj,
     const struct entry *restrict sort_i, const struct entry *restrict sort_j,
-    const float dx_max, const float rshift, float *max_di, float *max_dj,
-    int *init_pi, int *init_pj, const struct engine *e) {
+    const float dx_max, const float rshift, const double hi_max,
+    const double hj_max, const double di_max, const double dj_min,
+    float *max_di, float *max_dj, int *init_pi, int *init_pj,
+    const struct engine *e) {
 
-  struct part *restrict parts_i = ci->parts;
-  struct part *restrict parts_j = cj->parts;
-  struct part *p = &parts_i[sort_i[0].i];
-
-  float h, d;
-
-  /* Get the distance of the last pi and the first pj on the sorted axis.*/
-  const float di_max = sort_i[ci->count - 1].d - rshift;
-  const float dj_min = sort_j[0].d;
+  const struct part *restrict parts_i = ci->parts;
+  const struct part *restrict parts_j = cj->parts;
 
   int first_pi = 0, last_pj = cj->count - 1;
 
@@ -304,19 +303,19 @@ __attribute__((always_inline)) INLINE static void populate_max_d_no_cache(
   /* Populate max_di with distances. */
   int active_id = ci->count - 1;
   for (int k = ci->count - 1; k >= 0; k--) {
-    p = &parts_i[sort_i[k].i];
-    h = p->h;
-    d = sort_i[k].d + h * kernel_gamma + dx_max - rshift;
+    const struct part *pi = &parts_i[sort_i[k].i];
+    const float d = sort_i[k].d + dx_max;
 
-    max_di[k] = d;
+    // max_di[k] = d + h * kernel_gamma - rshift;
+    max_di[k] = d + hi_max;
 
     /* If the particle is out of range set the index to
      * the last active particle within range. */
-    if (d < dj_min) {
+    if (d + hi_max < dj_min) {
       first_pi = active_id;
       break;
     } else {
-      if (part_is_active(p, e)) active_id = k;
+      if (part_is_active(pi, e)) active_id = k;
     }
   }
 
@@ -329,19 +328,20 @@ __attribute__((always_inline)) INLINE static void populate_max_d_no_cache(
   /* Populate max_dj with distances. */
   active_id = 0;
   for (int k = 0; k < cj->count; k++) {
-    p = &parts_j[sort_j[k].i];
-    h = p->h;
-    d = sort_j[k].d - h * kernel_gamma - dx_max - rshift;
+    const struct part *pj = &parts_j[sort_j[k].i];
+    const float d = sort_j[k].d - dx_max;
 
-    max_dj[k] = d;
+    /*TODO: don't think rshift should be taken off here, waiting on Pedro. */
+    // max_dj[k] = d - h * kernel_gamma - rshift;
+    max_dj[k] = d - hj_max;
 
     /* If the particle is out of range set the index to
      * the last active particle within range. */
-    if (d > di_max) {
+    if (d - hj_max > di_max) {
       last_pj = active_id;
       break;
     } else {
-      if (part_is_active(p, e)) active_id = k;
+      if (part_is_active(pj, e)) active_id = k;
     }
   }
 
@@ -611,9 +611,12 @@ __attribute__((always_inline)) INLINE void runner_doself1_density_vec(
  * @param r The #runner.
  * @param ci The first #cell.
  * @param cj The second #cell.
+ * @param sid The direction of the pair
+ * @param shift The shift vector to apply to the particles in ci.
  */
 void runner_dopair1_density_vec(struct runner *r, struct cell *ci,
-                                struct cell *cj) {
+                                struct cell *cj, const int sid,
+                                const double *shift) {
 
 #ifdef WITH_VECTORIZATION
   const struct engine *restrict e = r->e;
@@ -622,22 +625,6 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci,
 
   TIMER_TIC;
 
-  /* Anything to do here? */
-  if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
-
-  if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e))
-    error("Interacting undrifted cells.");
-
-  /* Get the sort ID. */
-  double shift[3] = {0.0, 0.0, 0.0};
-  const int sid = space_getsid(e->s, &ci, &cj, shift);
-
-  /* Have the cells been sorted? */
-  if (!(ci->sorted & (1 << sid)) || ci->dx_max_sort > space_maxreldx * ci->dmin)
-    runner_do_sort(r, ci, (1 << sid), 1);
-  if (!(cj->sorted & (1 << sid)) || cj->dx_max_sort > space_maxreldx * cj->dmin)
-    runner_do_sort(r, cj, (1 << sid), 1);
-
   /* Get the cutoff shift. */
   double rshift = 0.0;
   for (int k = 0; k < 3; k++) rshift += shift[k] * runner_shift[sid][k];
@@ -726,8 +713,9 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci,
   /* Find particles maximum distance into cj, max_di[] and ci, max_dj[]. */
   /* Also find the first pi that interacts with any particle in cj and the last
    * pj that interacts with any particle in ci. */
-  populate_max_d_no_cache(ci, cj, sort_i, sort_j, dx_max, rshift, max_di,
-                          max_dj, &first_pi, &last_pj, e);
+  populate_max_d_no_cache(ci, cj, sort_i, sort_j, dx_max, rshift, hi_max,
+                          hj_max, di_max, dj_min, max_di, max_dj, &first_pi,
+                          &last_pj, e);
 
   /* Find the maximum index into cj that is required by a particle in ci. */
   /* Find the maximum index into ci that is required by a particle in cj. */
@@ -770,13 +758,22 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci,
 
   if (cell_is_active(ci, e)) {
 
-    /* Loop over the parts in ci. */
+    /* Loop over the parts in ci until nothing is within range in cj. */
     for (int pid = count_i - 1; pid >= first_pi_loop && max_ind_j >= 0; pid--) {
 
       /* Get a hold of the ith part in ci. */
       struct part *restrict pi = &parts_i[sort_i[pid].i];
       if (!part_is_active(pi, e)) continue;
 
+      /* Set the cache index. */
+      int ci_cache_idx = pid - first_pi_align;
+
+      /* Skip this particle if no particle in cj is within range of it. */
+      const float hi = ci_cache->h[ci_cache_idx];
+      const double di_test =
+          sort_i[pid].d + hi * kernel_gamma + dx_max - rshift;
+      if (di_test < dj_min) continue;
+
       /* Determine the exit iteration of the interaction loop. */
       dj = sort_j[max_ind_j].d;
       while (max_ind_j > 0 && max_di[pid] < dj) {
@@ -786,10 +783,6 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci,
       }
       int exit_iteration = max_ind_j + 1;
 
-      /* Set the cache index. */
-      int ci_cache_idx = pid - first_pi_align;
-
-      const float hi = ci_cache->h[ci_cache_idx];
       const float hig2 = hi * hi * kernel_gamma2;
 
       vector pix, piy, piz;
@@ -903,13 +896,24 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci,
   }
 
   if (cell_is_active(cj, e)) {
-    /* Loop over the parts in cj. */
+
+    /* Loop over the parts in cj until nothing is within range in ci. */
     for (int pjd = 0; pjd <= last_pj_loop && max_ind_i < count_i; pjd++) {
 
       /* Get a hold of the jth part in cj. */
       struct part *restrict pj = &parts_j[sort_j[pjd].i];
       if (!part_is_active(pj, e)) continue;
 
+      /* Set the cache index. */
+      int cj_cache_idx = pjd;
+
+      /*TODO: rshift term. */
+      /* Skip this particle if no particle in ci is within range of it. */
+      const float hj = cj_cache->h[cj_cache_idx];
+      const double dj_test =
+          sort_j[pjd].d - hj * kernel_gamma - dx_max - rshift;
+      if (dj_test > di_max) continue;
+
       /* Determine the exit iteration of the interaction loop. */
       di = sort_i[max_ind_i].d;
       while (max_ind_i < count_i - 1 && max_dj[pjd] > di) {
@@ -919,10 +923,6 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci,
       }
       int exit_iteration = max_ind_i;
 
-      /* Set the cache index. */
-      int cj_cache_idx = pjd;
-
-      const float hj = cj_cache->h[cj_cache_idx];
       const float hjg2 = hj * hj * kernel_gamma2;
 
       vector pjx, pjy, pjz;

src/runner_doiact_vec.h

@@ -35,8 +35,8 @@
 
 /* Function prototypes. */
 void runner_doself1_density_vec(struct runner *r, struct cell *restrict c);
-void runner_doself1_density_vec_2(struct runner *r, struct cell *restrict c);
 void runner_dopair1_density_vec(struct runner *r, struct cell *restrict ci,
-                                struct cell *restrict cj);
+                                struct cell *restrict cj, const int sid,
+                                const double *shift);
 
 #endif /* SWIFT_RUNNER_VEC_H */

tests/Makefile.am

@@ -94,6 +94,6 @@ EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
 	     test27cells.sh test27cellsPerturbed.sh testParser.sh \
 	     test125cells.sh test125cellsPerturbed.sh testParserInput.yaml difffloat.py \
 	     tolerance_125_normal.dat tolerance_125_perturbed.dat \
-             tolerance_27_normal.dat tolerance_27_perturbed.dat \
+             tolerance_27_normal.dat tolerance_27_perturbed.dat tolerance_27_perturbed_h.dat \
 	     tolerance_pair_normal.dat tolerance_pair_perturbed.dat \
 	     fft_params.yml

tests/difffloat.py

@@ -46,7 +46,13 @@ if len(sys.argv) == 6:
     ignoreSmallRhoDh = int(sys.argv[5])
 else:
     ignoreSmallRhoDh = 0
-    
+
+# Get the particle properties being compared from the header.
+with open(file1, 'r') as f:
+  line = f.readline()
+  if 'ID' in line:
+    part_props = line.split()[1:]
+
 data1 = loadtxt(file1)
 data2 = loadtxt(file2)
 if fileTol != "":
@@ -100,7 +106,7 @@ for i in range(n_lines_to_check):
             rel_diff = 0.
 
         if( abs_diff > 1.1*absTol[j]):
-            print "Absolute difference larger than tolerance (%e) for particle %d, column %d:"%(absTol[j], i,j)
+            print "Absolute difference larger than tolerance (%e) for particle %d, column %s:"%(absTol[j], data1[i,0], part_props[j])
             print "%10s:           a = %e"%("File 1", data1[i,j])
             print "%10s:           b = %e"%("File 2", data2[i,j])
             print "%10s:       |a-b| = %e"%("Difference", abs_diff)
@@ -113,7 +119,7 @@ for i in range(n_lines_to_check):
         if ignoreSmallRhoDh and j == 8 and abs(data1[i,j]) < 2e-4: continue
         
         if( rel_diff > 1.1*relTol[j]):
-            print "Relative difference larger than tolerance (%e) for particle %d, column %d:"%(relTol[j], i,j)
+            print "Relative difference larger than tolerance (%e) for particle %d, column %s:"%(relTol[j], data1[i,0], part_props[j])
             print "%10s:           a = %e"%("File 1", data1[i,j])
             print "%10s:           b = %e"%("File 2", data2[i,j])
             print "%10s: |a-b|/|a+b| = %e"%("Difference", rel_diff)

tests/test125cells.c

@@ -432,6 +432,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
 
 /* Just a forward declaration... */
 void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
+void runner_dopair1_branch_density(struct runner *r, struct cell *ci,
+                                   struct cell *cj);
 void runner_doself1_density(struct runner *r, struct cell *ci);
 void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
 void runner_doself2_force(struct runner *r, struct cell *ci);
@@ -621,6 +623,14 @@ int main(int argc, char *argv[]) {
 /* Do the density calculation */
 #if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
 
+    /* Initialise the particle cache. */
+#ifdef WITH_VECTORIZATION
+    runner.ci_cache.count = 0;
+    cache_init(&runner.ci_cache, 512);
+    runner.cj_cache.count = 0;
+    cache_init(&runner.cj_cache, 512);
+#endif
+
     /* Run all the pairs (only once !)*/
     for (int i = 0; i < 5; i++) {
       for (int j = 0; j < 5; j++) {
@@ -643,7 +653,7 @@ int main(int argc, char *argv[]) {
 
                 struct cell *cj = cells[iii * 25 + jjj * 5 + kkk];
 
-                if (cj > ci) runner_dopair1_density(&runner, ci, cj);
+                if (cj > ci) runner_dopair1_branch_density(&runner, ci, cj);
               }
             }
           }

tests/test27cells.c

@@ -30,11 +30,9 @@
 /* Local headers. */
 #include "swift.h"
 
-#define ACC_THRESHOLD 1e-5
-
 #if defined(WITH_VECTORIZATION)
 #define DOSELF1 runner_doself1_density_vec
-#define DOPAIR1 runner_dopair1_density_vec
+#define DOPAIR1 runner_dopair1_branch_density
 #define DOSELF1_NAME "runner_doself1_density_vec"
 #define DOPAIR1_NAME "runner_dopair1_density_vec"
 #endif
@@ -45,7 +43,7 @@
 #endif
 
 #ifndef DOPAIR1
-#define DOPAIR1 runner_dopair1_density
+#define DOPAIR1 runner_dopair1_branch_density
 #define DOPAIR1_NAME "runner_dopair1_density"
 #endif
 
@@ -64,18 +62,20 @@ enum velocity_types {
  * @param offset The position of the cell offset from (0,0,0).
  * @param size The cell size.
  * @param h The smoothing length of the particles in units of the inter-particle
- *separation.
+ * separation.
  * @param density The density of the fluid.
  * @param partId The running counter of IDs.
  * @param pert The perturbation to apply to the particles in the cell in units
- *of the inter-particle separation.
+ * of the inter-particle separation.
  * @param vel The type of velocity field (0, random, divergent, rotating)
+ * @param h_pert The perturbation to apply to the smoothing length.
  */
 struct cell *make_cell(size_t n, double *offset, double size, double h,
                        double density, long long *partId, double pert,
-                       enum velocity_types vel) {
+                       enum velocity_types vel, double h_pert) {
   const size_t count = n * n * n;
   const double volume = size * size * size;
+  float h_max = 0.f;
   struct cell *cell = malloc(sizeof(struct cell));
   bzero(cell, sizeof(struct cell));
 
@@ -121,7 +121,11 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
             part->v[2] = 0.f;
             break;
         }
-        part->h = size * h / (float)n;
+        if (h_pert)
+          part->h = size * h * random_uniform(1.f, 1.1f) / (float)n;
+        else
+          part->h = size * h / (float)n;
+        h_max = fmaxf(h_max, part->h);
         part->id = ++(*partId);
 
 #if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
@@ -156,7 +160,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
 
   /* Cell properties */
   cell->split = 0;
-  cell->h_max = h;
+  cell->h_max = h_max;
   cell->count = count;
   cell->dx_max_part = 0.;
   cell->dx_max_sort = 0.;
@@ -288,33 +292,11 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
   fclose(file);
 }
 
-/**
- * @brief Compares the vectorised result against
- * the serial result of the interaction.
- *
- * @param serial_parts Particle array that has been interacted serially
- * @param vec_parts Particle array to be interacted using vectors
- * @param count No. of particles that have been interacted
- * @param threshold Level of accuracy needed
- *
- * @return Non-zero value if difference found, 0 otherwise
- */
-int check_results(struct part *serial_parts, struct part *vec_parts, int count,
-                  double threshold) {
-  int result = 0;
-
-  for (int i = 0; i < count; i++)
-    result += compare_particles(serial_parts[i], vec_parts[i], threshold);
-
-  return result;
-}
-
 /* Just a forward declaration... */
 void runner_doself1_density(struct runner *r, struct cell *ci);
 void runner_doself1_density_vec(struct runner *r, struct cell *ci);
-void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
-void runner_dopair1_density_vec(struct runner *r, struct cell *ci,
-                                struct cell *cj);
+void runner_dopair1_branch_density(struct runner *r, struct cell *ci,
+                                   struct cell *cj);
 
 /* And go... */
 int main(int argc, char *argv[]) {
@@ -322,8 +304,7 @@ int main(int argc, char *argv[]) {
   engine_pin();
   size_t runs = 0, particles = 0;
   double h = 1.23485, size = 1., rho = 1.;
-  double perturbation = 0.;
-  double threshold = ACC_THRESHOLD;
+  double perturbation = 0., h_pert = 0.;
   char outputFileNameExtension[200] = "";
   char outputFileName[200] = "";
   enum velocity_types vel = velocity_zero;
@@ -339,11 +320,14 @@ int main(int argc, char *argv[]) {
   srand(0);
 
   char c;
-  while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:a:")) != -1) {
+  while ((c = getopt(argc, argv, "m:s:h:p:n:r:t:d:f:v:")) != -1) {
     switch (c) {
       case 'h':
         sscanf(optarg, "%lf", &h);
         break;
+      case 'p':
+        sscanf(optarg, "%lf", &h_pert);
+        break;
       case 's':
         sscanf(optarg, "%lf", &size);
         break;
@@ -365,9 +349,6 @@ int main(int argc, char *argv[]) {
       case 'v':
         sscanf(optarg, "%d", (int *)&vel);
         break;
-      case 'a':
-        sscanf(optarg, "%lf", &threshold);
-        break;
       case '?':
         error("Unknown option.");
         break;
@@ -382,6 +363,7 @@ int main(int argc, char *argv[]) {
         "runner_doself1_density()."
         "\n\nOptions:"
         "\n-h DISTANCE=1.2348 - Smoothing length in units of <x>"
+        "\n-p                 - Random fractional change in h, h=h*random(1,p)"
         "\n-m rho             - Physical density in the cell"
         "\n-s size            - Physical size of the cell"
         "\n-d pert            - Perturbation to apply to the particles [0,1["
@@ -435,8 +417,9 @@ int main(int argc, char *argv[]) {
     for (int j = 0; j < 3; ++j) {
       for (int k = 0; k < 3; ++k) {
         double offset[3] = {i * size, j * size, k * size};
-        cells[i * 9 + j * 3 + k] = make_cell(particles, offset, size, h, rho,
-                                             &partId, perturbation, vel);
+        cells[i * 9 + j * 3 + k] =
+            make_cell(particles, offset, size, h, rho, &partId, perturbation,
+                      vel, h_pert);
 
         runner_do_drift_part(&runner, cells[i * 9 + j * 3 + k], 0);
 
@@ -504,10 +487,6 @@ int main(int argc, char *argv[]) {
     }
   }
 
-  /* Store the vectorised particle results. */
-  struct part vec_parts[main_cell->count];
-  for (int i = 0; i < main_cell->count; i++) vec_parts[i] = main_cell->parts[i];
-
   /* Output timing */
   ticks corner_time = timings[0] + timings[2] + timings[6] + timings[8] +
                       timings[18] + timings[20] + timings[24] + timings[26];
@@ -552,10 +531,6 @@ int main(int argc, char *argv[]) {