Initial setup of a new test case with 125 cells.

.gitignore

@@ -39,6 +39,7 @@ tests/swift_dopair_standard.dat
 tests/brute_force_perturbed.dat
 tests/swift_dopair_perturbed.dat
 tests/test27cells
+tests/test125cells
 tests/brute_force_27_standard.dat
 tests/swift_dopair_27_standard.dat
 tests/brute_force_27_perturbed.dat
@@ -54,6 +55,7 @@ tests/testParser
 tests/parser_output.yml
 tests/test27cells.sh
 tests/test27cellsPerturbed.sh
+tests/test125cells.sh
 tests/testPair.sh
 tests/testPairPerturbed.sh
 tests/testParser.sh

configure.ac

@@ -476,6 +476,7 @@ AC_CONFIG_FILES([tests/testPair.sh], [chmod +x tests/testPair.sh])
 AC_CONFIG_FILES([tests/testPairPerturbed.sh], [chmod +x tests/testPairPerturbed.sh])
 AC_CONFIG_FILES([tests/test27cells.sh], [chmod +x tests/test27cells.sh])
 AC_CONFIG_FILES([tests/test27cellsPerturbed.sh], [chmod +x tests/test27cellsPerturbed.sh])
+AC_CONFIG_FILES([tests/test125cells.sh], [chmod +x tests/test125cells.sh])
 AC_CONFIG_FILES([tests/testParser.sh], [chmod +x tests/testParser.sh])
 
 # Report general configuration.

tests/Makefile.am

@@ -22,11 +22,13 @@ AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
 
 # List of programs and scripts to run in the test suite
 TESTS = testGreetings testReading.sh testSingle testPair.sh testPairPerturbed.sh \
-	test27cells.sh test27cellsPerturbed.sh testParser.sh testKernel testSPHStep
+	test27cells.sh test27cellsPerturbed.sh testParser.sh testKernel testSPHStep \
+	test125cells.sh
 
 # List of test programs to compile
 check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
-		 testSPHStep testPair test27cells testParser testKernel testNonSymInt testSymInt
+		 testSPHStep testPair test27cells test125cells testParser testKernel \ 
+                 testNonSymInt testSymInt
 
 # Sources for the individual programs
 testGreetings_SOURCES = testGreetings.c
@@ -43,6 +45,8 @@ testPair_SOURCES = testPair.c
 
 test27cells_SOURCES = test27cells.c
 
+test125cells_SOURCES = test125cells.c
+
 testParser_SOURCES = testParser.c
 
 testKernel_SOURCES = testKernel.c
@@ -54,4 +58,4 @@ testSymInt_SOURCES = testSymInt.c
 # Files necessary for distribution
 EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
 	     test27cells.sh test27cellsPerturbed.sh tolerance.dat testParser.sh \
-	     testParserInput.yaml
+	     test125cells.sh testParserInput.yaml

tests/test125cells.c

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include <fenv.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+#include "swift.h"
+
+enum velocity_types {
+  velocity_zero,
+  velocity_random,
+  velocity_divergent,
+  velocity_rotating
+};
+
+/**
+ * @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 vel The type of velocity field (0, random, divergent, rotating)
+ */
+struct cell *make_cell(size_t n, double *offset, double size, double h,
+                       double density, long long *partId, double pert,
+                       enum velocity_types vel) {
+  const size_t count = n * n * n;
+  const double volume = size * size * size;
+  struct cell *cell = malloc(sizeof(struct cell));
+  bzero(cell, sizeof(struct cell));
+
+  if (posix_memalign((void **)&cell->parts, part_align,
+                     count * sizeof(struct part)) != 0) {
+    error("couldn't allocate particles, no. of particles: %d", (int)count);
+  }
+  bzero(cell->parts, count * sizeof(struct part));
+
+  /* Construct the parts */
+  struct part *part = cell->parts;
+  for (size_t x = 0; x < n; ++x) {
+    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;
+        part->x[1] =
+            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;
+        switch (vel) {
+          case velocity_zero:
+            part->v[0] = 0.f;
+            part->v[1] = 0.f;
+            part->v[2] = 0.f;
+            break;
+          case velocity_random:
+            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);
+            break;
+          case velocity_divergent:
+            part->v[0] = part->x[0] - 1.5 * size;
+            part->v[1] = part->x[1] - 1.5 * size;
+            part->v[2] = part->x[2] - 1.5 * size;
+            break;
+          case velocity_rotating:
+            part->v[0] = part->x[1];
+            part->v[1] = -part->x[0];
+            part->v[2] = 0.f;
+            break;
+        }
+        part->h = size * h / (float)n;
+        part->id = ++(*partId);
+        part->mass = density * volume / count;
+        part->ti_begin = 0;
+        part->ti_end = 1;
+        ++part;
+      }
+    }
+  }
+
+  /* Cell properties */
+  cell->split = 0;
+  cell->h_max = h;
+  cell->count = count;
+  cell->dx_max = 0.;
+  cell->h[0] = size;
+  cell->h[1] = size;
+  cell->h[2] = size;
+  cell->loc[0] = offset[0];
+  cell->loc[1] = offset[1];
+  cell->loc[2] = offset[2];
+
+  cell->ti_end_min = 1;
+  cell->ti_end_max = 1;
+
+  shuffle_particles(cell->parts, cell->count);
+
+  cell->sorted = 0;
+  cell->sort = NULL;
+  cell->sortsize = 0;
+
+  return cell;
+}
+
+void clean_up(struct cell *ci) {
+  free(ci->parts);
+  free(ci->sort);
+  free(ci);
+}
+
+/**
+ * @brief Initializes all particles field to be ready for a density calculation
+ */
+void zero_particle_fields(struct cell *c) {
+  for (size_t pid = 0; pid < c->count; pid++) {
+    c->parts[pid].rho = 0.f;
+    c->parts[pid].rho_dh = 0.f;
+    hydro_init_part(&c->parts[pid]);
+  }
+}
+
+/**
+ * @brief Ends the loop by adding the appropriate coefficients
+ */
+void end_calculation(struct cell *c) {
+  for (size_t pid = 0; pid < c->count; pid++) {
+    hydro_end_density(&c->parts[pid], 1);
+  }
+}
+
+/**
+ * @brief Dump all the particles to a file
+ */
+void dump_particle_fields(char *fileName, struct cell *main_cell,
+                          struct cell **cells) {
+  FILE *file = fopen(fileName, "w");
+
+  /* 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");
+
+  fprintf(file, "# Main cell --------------------------------------------\n");
+
+  /* Write main cell */
+  for (size_t pid = 0; pid < main_cell->count; pid++) {
+    fprintf(file,
+            "%6llu %10f %10f %10f %10f %10f %10f %13e %13e %13e %13e %13e "
+            "%13e %13e %13e\n",
+            main_cell->parts[pid].id, main_cell->parts[pid].x[0],
+            main_cell->parts[pid].x[1], main_cell->parts[pid].x[2],
+            main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
+            main_cell->parts[pid].v[2], main_cell->parts[pid].rho,
+            main_cell->parts[pid].rho_dh, main_cell->parts[pid].density.wcount,
+            main_cell->parts[pid].density.wcount_dh,
+#if defined(GADGET2_SPH)
+            main_cell->parts[pid].div_v, main_cell->parts[pid].density.rot_v[0],
+            main_cell->parts[pid].density.rot_v[1],
+            main_cell->parts[pid].density.rot_v[2]
+#elif defined(DEFAULT_SPH)
+            main_cell->parts[pid].density.div_v,
+            main_cell->parts[pid].density.rot_v[0],
+            main_cell->parts[pid].density.rot_v[1],
+            main_cell->parts[pid].density.rot_v[2]
+#else
+            0., 0., 0., 0.
+#endif
+            );
+  }
+
+  /* Write all other cells */
+  for (int i = 0; i < 3; ++i) {
+    for (int j = 0; j < 3; ++j) {
+      for (int k = 0; k < 3; ++k) {
+        struct cell *cj = cells[i * 9 + j * 3 + k];
+        if (cj == main_cell) continue;
+
+        fprintf(file,
+                "# Offset: [%2d %2d %2d] -----------------------------------\n",
+                i - 1, j - 1, k - 1);
+
+        for (size_t 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], cj->parts[pjd].rho, cj->parts[pjd].rho_dh,
+              cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
+#if defined(GADGET2_SPH)
+              cj->parts[pjd].div_v, cj->parts[pjd].density.rot_v[0],
+              cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
+#elif defined(DEFAULT_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
+              );
+        }
+      }
+    }
+  }
+  fclose(file);
+}
+
+/* Just a forward declaration... */
+void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
+void runner_doself1_density(struct runner *r, struct cell *ci);
+
+/* And go... */
+int main(int argc, char *argv[]) {
+  size_t runs = 0, particles = 0;
+  double h = 1.2348, size = 1., rho = 1.;
+  double perturbation = 0.;
+  char outputFileNameExtension[200] = "";
+  char outputFileName[200] = "";
+  int vel = velocity_zero;
+
+  /* Initialize CPU frequency, this also starts time. */
+  unsigned long long cpufreq = 0;
+  clocks_set_cpufreq(cpufreq);
+
+  /* Get some randomness going */
+  srand(0);
+
+  char c;
+  while ((c = getopt(argc, argv, "m:s:h:p:r:t:d:f:v:")) != -1) {
+    switch (c) {
+      case 'h':
+        sscanf(optarg, "%lf", &h);
+        break;
+      case 's':
+        sscanf(optarg, "%lf", &size);
+        break;
+      case 'p':
+        sscanf(optarg, "%zu", &particles);
+        break;
+      case 'r':
+        sscanf(optarg, "%zu", &runs);
+        break;
+      case 'd':
+        sscanf(optarg, "%lf", &perturbation);
+        break;
+      case 'm':
+        sscanf(optarg, "%lf", &rho);
+        break;
+      case 'f':
+        strcpy(outputFileNameExtension, optarg);
+        break;
+      case 'v':
+        sscanf(optarg, "%d", &vel);
+        break;
+      case '?':
+        error("Unknown option.");
+        break;
+    }
+  }
+
+  if (h < 0 || particles == 0 || runs == 0) {
+    printf(
+        "\nUsage: %s -p 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-h DISTANCE=1.2348 - Smoothing length in units of <x>"
+        "\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["
+        "\n-v type (0,1,2,3)  - Velocity field: (zero, random, divergent, "
+        "rotating)"
+        "\n-f fileName        - Part of the file name used to save the dumps\n",
+        argv[0]);
+    exit(1);
+  }
+
+  /* Help users... */
+  message("Smoothing length: h = %f", h * size);
+  message("Kernel:               %s", kernel_name);
+  message("Neighbour target: N = %f",
+          h * h * h * 4.0 * M_PI * kernel_gamma3 / 3.0);
+  message("Density target: rho = %f", rho);
+  message("div_v target:   div = %f", vel == 2 ? 3.f : 0.f);
+  message("curl_v target: curl = [0., 0., %f]", vel == 3 ? -2.f : 0.f);
+  printf("\n");
+
+  /* Build the infrastructure */
+  struct space space;
+  space.periodic = 0;
+  space.h_max = h;
+
+  struct engine engine;
+  engine.s = &space;
+  engine.time = 0.1f;
+  engine.ti_current = 1;
+
+  struct runner runner;
+  runner.e = &engine;
+
+  /* Construct some cells */
+  struct cell *cells[27];
+  struct cell *main_cell;
+  static long long partId = 0;
+  for (int i = 0; i < 3; ++i) {
+    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);
+
+        runner_do_sort(&runner, cells[i * 9 + j * 3 + k], 0x1FFF, 0);
+      }
+    }
+  }
+
+  /* Store the main cell for future use */
+  main_cell = cells[13];
+
+  ticks time = 0;
+  for (size_t i = 0; i < runs; ++i) {
+    /* Zero the fields */
+    for (int j = 0; j < 27; ++j) zero_particle_fields(cells[j]);
+
+    const ticks tic = getticks();
+
+#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
+
+    /* Run all the pairs */
+    for (int j = 0; j < 27; ++j)
+      if (cells[j] != main_cell)
+        runner_dopair1_density(&runner, main_cell, cells[j]);
+
+    /* And now the self-interaction */
+    runner_doself1_density(&runner, main_cell);
+
+#endif
+
+    const ticks toc = getticks();
+    time += toc - tic;
+
+    /* Let's get physical ! */
+    end_calculation(main_cell);
+
+    /* Dump if necessary */
+    if (i % 50 == 0) {
+      sprintf(outputFileName, "swift_dopair_27_%s.dat",
+              outputFileNameExtension);
+      dump_particle_fields(outputFileName, main_cell, cells);
+    }
+  }
+
+  /* Output timing */
+  message("SWIFT calculation took       : %15lli ticks.", time / runs);
+
+  /* Now perform a brute-force version for accuracy tests */
+
+  /* Zero the fields */
+  for (int i = 0; i < 27; ++i) zero_particle_fields(cells[i]);
+
+  const ticks tic = getticks();
+
+#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
+
+  /* Run all the brute-force pairs */
+  for (int j = 0; j < 27; ++j)
+    if (cells[j] != main_cell) pairs_all_density(&runner, main_cell, cells[j]);
+
+  /* And now the self-interaction */
+  self_all_density(&runner, main_cell);
+
+#endif
+
+  const ticks toc = getticks();
+
+  /* Let's get physical ! */
+  end_calculation(main_cell);
+
+  /* Dump */
+  sprintf(outputFileName, "brute_force_27_%s.dat", outputFileNameExtension);
+  dump_particle_fields(outputFileName, main_cell, cells);
+
+  /* Output timing */
+  message("Brute force calculation took : %15lli ticks.", toc - tic);
+
+  /* Clean things to make the sanitizer happy ... */
+  for (int i = 0; i < 27; ++i) clean_up(cells[i]);
+
+  return 0;
+}

tests/test125cells.sh.in

+#!/bin/bash
+rm brute_force_125_standard.dat swift_dopair_125_standard.dat
+
+./test125cells -p 4 -r 1 -d 0 -f standard
+
+#python @srcdir@/difffloat.py brute_force_125_standard.dat swift_dopair_125_standard.dat @srcdir@/tolerance.dat 4
+
+exit $?