Only tests cells with particles that are all active.

tests/testActivePair.c

@@ -28,14 +28,27 @@
 /* Local headers. */
 #include "swift.h"
 
-/* n is both particles per axis and box size:
- * particles are generated on a mesh with unit spacing
+/**
+ * @brief Constructs a cell and all of its particle in a valid state prior to
+ * a DOPAIR or DOSELF calcuation.
+ *
+ * @param n The cube root of the number of particles.
+ * @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.
+ * @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.
+ * @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, unsigned long long *partId,
-                       double pert) {
+                       double density, long long *partId, double pert,
+                       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));
 
@@ -51,26 +64,44 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
     for (size_t y = 0; y < n; ++y) {
       for (size_t z = 0; z < n; ++z) {
         part->x[0] =
-            offset[0] +
-            size * (x + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
+          offset[0] +
+          size * (x + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
         part->x[1] =
-            offset[1] +
-            size * (y + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
+          offset[1] +
+          size * (y + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
         part->x[2] =
-            offset[2] +
-            size * (z + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
+          offset[2] +
+          size * (z + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
         part->v[0] = random_uniform(-0.05, 0.05);
         part->v[1] = random_uniform(-0.05, 0.05);
         part->v[2] = random_uniform(-0.05, 0.05);
-        part->h = size * h / (float)n;
+
+        if (h_pert)
+          part->h = size * h * random_uniform(1.f, h_pert) / (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)
         part->conserved.mass = density * volume / count;
+
+#ifdef SHADOWFAX_SPH
+        double anchor[3] = {0., 0., 0.};
+        double side[3] = {1., 1., 1.};
+        voronoi_cell_init(&part->cell, part->x, anchor, side);
+#endif
+
 #else
         part->mass = density * volume / count;
 #endif
+
+#if defined(HOPKINS_PE_SPH)
+        part->entropy = 1.f;
+        part->entropy_one_over_gamma = 1.f;
+#endif
+
         part->time_bin = 1;
-        //part->time_bin = num_time_bins + 1;
 
 #ifdef SWIFT_DEBUG_CHECKS
         part->ti_drift = 8;
@@ -84,13 +115,13 @@ 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.;
-  cell->width[0] = n;
-  cell->width[1] = n;
-  cell->width[2] = n;
+  cell->width[0] = size;
+  cell->width[1] = size;
+  cell->width[2] = size;
   cell->loc[0] = offset[0];
   cell->loc[1] = offset[1];
   cell->loc[2] = offset[2];
@@ -98,6 +129,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
   cell->ti_old_part = 8;
   cell->ti_end_min = 8;
   cell->ti_end_max = 8;
+  cell->ti_sort = 8;
 
   shuffle_particles(cell->parts, cell->count);
 
@@ -122,6 +154,19 @@ void zero_particle_fields(struct cell *c) {
   }
 }
 
+/**
+ * @brief Ends the loop by adding the appropriate coefficients
+ */
+void end_calculation(struct cell *c) {
+  for (int pid = 0; pid < c->count; pid++) {
+    hydro_end_density(&c->parts[pid]);
+
+    /* Recover the common "Neighbour number" definition */
+    c->parts[pid].density.wcount *= pow_dimension(c->parts[pid].h);
+    c->parts[pid].density.wcount *= kernel_norm;
+  }
+}
+
 /**
  * @brief Dump all the particles to a file
  */
@@ -130,57 +175,24 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
 
   /* Write header */
   fprintf(file,
-          "# %4s %10s %10s %10s %10s %10s %10s %13s %13s %13s %13s %13s "
-          "%13s %13s %13s\n",
-          "ID", "pos_x", "pos_y", "pos_z", "v_x", "v_y", "v_z", "rho", "rho_dh",
-          "wcount", "wcount_dh", "div_v", "curl_vx", "curl_vy", "curl_vz");
+          "# %4s %13s\n", "ID", "wcount");
 
   fprintf(file, "# ci --------------------------------------------\n");
 
   for (int pid = 0; pid < ci->count; pid++) {
     fprintf(file,
-            "%6llu %10f %10f %10f %10f %10f %10f %13e %13e %13e %13e %13e "
-            "%13e %13e %13e\n",
-            ci->parts[pid].id, ci->parts[pid].x[0], ci->parts[pid].x[1],
-            ci->parts[pid].x[2], ci->parts[pid].v[0], ci->parts[pid].v[1],
-            ci->parts[pid].v[2], hydro_get_density(&ci->parts[pid]),
-#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
-            0.f,
-#else
-            ci->parts[pid].density.rho_dh,
-#endif
-            ci->parts[pid].density.wcount, ci->parts[pid].density.wcount_dh,
-#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
-            ci->parts[pid].density.div_v, ci->parts[pid].density.rot_v[0],
-            ci->parts[pid].density.rot_v[1], ci->parts[pid].density.rot_v[2]
-#else
-            0., 0., 0., 0.
-#endif
-            );
+            "%6llu %13e\n",
+            ci->parts[pid].id,
+            ci->parts[pid].density.wcount);
   }
 
   fprintf(file, "# cj --------------------------------------------\n");
 
   for (int pjd = 0; pjd < cj->count; pjd++) {
     fprintf(file,
-            "%6llu %10f %10f %10f %10f %10f %10f %13e %13e %13e %13e %13e "
-            "%13e %13e %13e\n",
-            cj->parts[pjd].id, cj->parts[pjd].x[0], cj->parts[pjd].x[1],
-            cj->parts[pjd].x[2], cj->parts[pjd].v[0], cj->parts[pjd].v[1],
-            cj->parts[pjd].v[2], hydro_get_density(&cj->parts[pjd]),
-#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
-            0.f,
-#else
-            cj->parts[pjd].density.rho_dh,
-#endif
-            cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
-#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
-            cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0],
-            cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
-#else
-            0., 0., 0., 0.
-#endif
-            );
+            "%6llu %13e\n",
+            cj->parts[pjd].id,
+            cj->parts[pjd].density.wcount);
   }
 
   fclose(file);
@@ -192,19 +204,64 @@ void runner_doself1_density_vec(struct runner *r, struct cell *ci);
 void runner_dopair1_branch_density(struct runner *r, struct cell *ci,
                                    struct cell *cj);
 
+void test_pair_interactions(struct runner runner,
+                              struct cell **ci, struct cell **cj, char *swiftOutputFileName,
+                              char *bruteForceOutputFileName) {
+
+  runner_do_sort(&runner, *ci, 0x1FFF, 0, 0);
+  runner_do_sort(&runner, *cj, 0x1FFF, 0, 0);
+
+  /* Zero the fields */
+  zero_particle_fields(*ci);
+  zero_particle_fields(*cj);
+
+#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 the test */
+  runner_dopair1_branch_density(&runner, *ci, *cj);
+
+  /* Let's get physical ! */
+  end_calculation(*ci);
+  end_calculation(*cj);
+  
+  /* Dump if necessary */
+  dump_particle_fields(swiftOutputFileName, *ci, *cj);
+
+  /* Now perform a brute-force version for accuracy tests */
+
+  /* Zero the fields */
+  zero_particle_fields(*ci);
+  zero_particle_fields(*cj);
+
+  /* Run the brute-force test */
+  pairs_all_density(&runner, *ci, *cj);
+
+  /* Let's get physical ! */
+  end_calculation(*ci);
+  end_calculation(*cj);
+
+  dump_particle_fields(bruteForceOutputFileName, *ci, *cj);
+
+}
+
 int main(int argc, char *argv[]) {
-  size_t particles = 0, runs = 0, volume, type = 0;
-  double offset[3] = {0, 0, 0}, h = 1.1255, size = 1., rho = 1.;
-  double perturbation = 0.1;
+  size_t particles = 0, runs = 0, type = 0;
+  double offset[3] = {0, 0, 0}, h = 1.23485, size = 1., rho = 1.;
+  double perturbation = 0.1, h_pert = 1.1;
   struct cell *ci, *cj;
   struct space space;
   struct engine engine;
   struct runner runner;
   char c;
-  static unsigned long long partId = 0;
+  static long long partId = 0;
   char outputFileNameExtension[200] = "";
-  char outputFileName[200] = "";
-  ticks tic, toc, time;
+  char swiftOutputFileName[200] = "";
+  char bruteForceOutputFileName[200] = "";
 
   /* Initialize CPU frequency, this also starts time. */
   unsigned long long cpufreq = 0;
@@ -215,12 +272,15 @@ int main(int argc, char *argv[]) {
 
   srand(0);
 
-  while ((c = getopt(argc, argv, "h:p:r:t:d:f:")) != -1) {
+  while ((c = getopt(argc, argv, "h:p:n:r:t:d:f:")) != -1) {
     switch (c) {
       case 'h':
         sscanf(optarg, "%lf", &h);
         break;
       case 'p':
+        sscanf(optarg, "%lf", &h_pert);
+        break;
+      case 'n':
         sscanf(optarg, "%zu", &particles);
         break;
       case 'r':
@@ -243,12 +303,13 @@ int main(int argc, char *argv[]) {
 
   if (h < 0 || particles == 0 || runs == 0 || type > 2) {
     printf(
-        "\nUsage: %s -p PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
+        "\nUsage: %s -n PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
         "\nGenerates a cell pair, filled with particles on a Cartesian grid."
         "\nThese are then interacted using runner_dopair1_density."
         "\n\nOptions:"
         "\n-t TYPE=0          - cells share face (0), edge (1) or corner (2)"
-        "\n-h DISTANCE=1.1255 - smoothing length"
+        "\n-h DISTANCE=1.2348 - smoothing length"
+        "\n-p                 - Random fractional change in h, h=h*random(1,p)"
         "\n-d pert            - perturbation to apply to the particles [0,1["
         "\n-f fileName        - part of the file name used to save the dumps\n",
         argv[0]);
@@ -263,62 +324,17 @@ int main(int argc, char *argv[]) {
   engine.max_active_bin = num_time_bins;
   runner.e = &engine;
 
-  volume = particles * particles * particles;
-  message("particles: %zu B\npositions: 0 B", 2 * volume * sizeof(struct part));
-
-  ci = make_cell(particles, offset, size, h, rho, &partId, perturbation);
+  /* Create output file names. */
+  sprintf(swiftOutputFileName, "swift_dopair_%s.dat",
+          outputFileNameExtension);
+  sprintf(bruteForceOutputFileName, "brute_force_%s.dat",
+          outputFileNameExtension);
+  
+  ci = make_cell(particles, offset, size, h, rho, &partId, perturbation, h_pert);
   for (size_t i = 0; i < type + 1; ++i) offset[i] = 1.;
-  cj = make_cell(particles, offset, size, h, rho, &partId, perturbation);
-
-  runner_do_sort(&runner, ci, 0x1FFF, 0, 0);
-  runner_do_sort(&runner, cj, 0x1FFF, 0, 0);
-
-  time = 0;
-  /* Zero the fields */
-  zero_particle_fields(ci);
-  zero_particle_fields(cj);
-
-#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
-
-  tic = getticks();
-
-  /* Run the test */
-  runner_dopair1_branch_density(&runner, ci, cj);
-
-  toc = getticks();
-  time += toc - tic;
-
-  /* Dump if necessary */
-  sprintf(outputFileName, "swift_dopair_%s.dat", outputFileNameExtension);
-  dump_particle_fields(outputFileName, ci, cj);
-
-  /* Output timing */
-  message("SWIFT calculation took       %lli ticks.", time / runs);
-
-  /* Now perform a brute-force version for accuracy tests */
-
-  /* Zero the fields */
-  zero_particle_fields(ci);
-  zero_particle_fields(cj);
-
-  tic = getticks();
-
-  /* Run the brute-force test */
-  pairs_all_density(&runner, ci, cj);
-
-  toc = getticks();
-
-  /* Dump */
-  sprintf(outputFileName, "brute_force_%s.dat", outputFileNameExtension);
-  dump_particle_fields(outputFileName, ci, cj);
+  cj = make_cell(particles, offset, size, h, rho, &partId, perturbation, h_pert);
 
-  /* Output timing */
-  message("Brute force calculation took %lli ticks.", toc - tic);
+  test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName);
 
   /* Clean things to make the sanitizer happy ... */
   clean_up(ci);

tests/testActivePair.sh.in

@@ -4,7 +4,7 @@ echo ""
 
 rm -f brute_force_pair_active.dat swift_dopair_active.dat
 
-./testActivePair -p 6 -r 1 -d 0 -f active
+./testActivePair -n 6 -r 1 -d 0 -f active
 
 python @srcdir@/difffloat.py brute_force_active.dat swift_dopair_active.dat @srcdir@/tolerance_pair_active.dat