Added a test of the FFT task#

src/runner.h

@@ -70,6 +70,7 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer);
 void runner_do_init(struct runner *r, struct cell *c, int timer);
 void runner_do_cooling(struct runner *r, struct cell *c, int timer);
 void runner_do_grav_external(struct runner *r, struct cell *c, int timer);
+void runner_do_grav_fft(struct runner *r, int timer);
 void *runner_main(void *data);
 void runner_do_unskip_mapper(void *map_data, int num_elements,
                              void *extra_data);

src/space.c

@@ -895,8 +895,9 @@ void space_rebuild(struct space *s, int verbose) {
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Verify that the links are correct */
-  part_verify_links(s->parts, s->gparts, s->sparts, nr_parts, nr_gparts,
-                    nr_sparts, verbose);
+  if ((nr_gparts > 0 && nr_parts > 0) || (nr_gparts > 0 && nr_sparts > 0))
+    part_verify_links(s->parts, s->gparts, s->sparts, nr_parts, nr_gparts,
+                      nr_sparts, verbose);
 #endif
 
   /* Hook the cells up to the parts. */
@@ -2952,47 +2953,3 @@ void space_clean(struct space *s) {
   free(s->gparts);
   free(s->sparts);
 }
-
-void space_dump_multipoles(struct space *s) {
-
-  hid_t h_file =
-      H5Fcreate("multipoles.hdf5", H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
-  hid_t h_grp =
-      H5Gcreate(h_file, "/Multipoles", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-
-  const hid_t h_space = H5Screate(H5S_SIMPLE);
-  const hid_t h_space2 = H5Screate(H5S_SIMPLE);
-  hsize_t shape[2] = {s->nr_cells, 3};
-  hsize_t shape2[1] = {s->nr_cells};
-  H5Sset_extent_simple(h_space, 2, shape, NULL);
-  H5Sset_extent_simple(h_space2, 1, shape2, NULL);
-
-  const hid_t h_data =
-      H5Dcreate(h_grp, "Coordinates", H5T_NATIVE_DOUBLE, h_space, H5P_DEFAULT,
-                H5P_DEFAULT, H5P_DEFAULT);
-
-  const hid_t h_data2 =
-      H5Dcreate(h_grp, "Potential", H5T_NATIVE_DOUBLE, h_space2, H5P_DEFAULT,
-                H5P_DEFAULT, H5P_DEFAULT);
-
-  double *temp = malloc(s->nr_cells * 3 * sizeof(double));
-  double *temp2 = malloc(s->nr_cells * 1 * sizeof(double));
-  for (int i = 0; i < s->nr_cells; ++i) {
-    temp[i * 3 + 0] = s->multipoles_top[i].CoM[0];
-    temp[i * 3 + 1] = s->multipoles_top[i].CoM[1];
-    temp[i * 3 + 2] = s->multipoles_top[i].CoM[2];
-    temp2[i] = s->multipoles_top[i].pot.F_000;
-  }
-
-  H5Dwrite(h_data, H5T_NATIVE_DOUBLE, h_space, H5S_ALL, H5P_DEFAULT, temp);
-  H5Dwrite(h_data2, H5T_NATIVE_DOUBLE, h_space2, H5S_ALL, H5P_DEFAULT, temp2);
-
-  free(temp);
-  free(temp2);
-  H5Dclose(h_data);
-  H5Dclose(h_data2);
-  H5Sclose(h_space);
-  H5Sclose(h_space2);
-  H5Gclose(h_grp);
-  H5Fclose(h_file);
-}

src/space.h

@@ -224,6 +224,5 @@ void space_check_timesteps(struct space *s);
 void space_replicate(struct space *s, int replicate, int verbose);
 void space_reset_task_counters(struct space *s);
 void space_clean(struct space *s);
-void space_dump_multipoles(struct space *s);
 
 #endif /* SWIFT_SPACE_H */

tests/Makefile.am

@@ -94,4 +94,5 @@ EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
 	     test27cells.sh test27cellsPerturbed.sh testParser.sh \
 	     test125cells.sh testParserInput.yaml difffloat.py \
 	     tolerance_125.dat tolerance_27_normal.dat tolerance_27_perturbed.dat \
-	     tolerance_pair_normal.dat tolerance_pair_perturbed.dat
+	     tolerance_pair_normal.dat tolerance_pair_perturbed.dat \
+	     fft_params.yml

tests/fft_params.yml

+Scheduler:
+  max_top_level_cells: 64
+  
+# 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.00001  # Softening length (in internal units).
+  a_smooth:              0.
+  r_cut:                 0.

tests/testFFT.c

@@ -18,8 +18,8 @@
  ******************************************************************************/
 
 /* Some standard headers. */
-#include <stdlib.h>
-#include <string.h>
+
+#include "../config.h"
 
 #ifndef HAVE_FFTW
 
@@ -27,19 +27,15 @@ int main() { return 0; }
 
 #else
 
-#include <fftw3.h>
+/* Some standard headers. */
+#include <fenv.h>
+#include <stdlib.h>
+#include <string.h>
 
 /* Includes. */
 #include "swift.h"
 
-const double G = 1.;
-
-const size_t N = 16;
-const size_t PMGRID = 8;
-
-// const double asmth = 2. * M_PI * const_gravity_a_smooth / boxSize;
-// const double asmth2 = asmth * asmth;
-// const double fact = G / (M_PI * boxSize) * (1. / (2. * boxSize / PMGRID));
+#define NEAREST(x) (((x) > 0.5) ? ((x)-1.) : (((x) < -0.5) ? ((x) + 1.) : (x)))
 
 int main() {
 
@@ -47,149 +43,77 @@ int main() {
   unsigned long long cpufreq = 0;
   clocks_set_cpufreq(cpufreq);
 
-  /* Simulation properties */
-  const size_t count = N * N * N;
-  const double boxSize = 1.;
-
-  /* Create some particles */
-  struct gpart* gparts = malloc(count * sizeof(struct gpart));
-  bzero(gparts, count * sizeof(struct gpart));
-  for (size_t i = 0; i < N; ++i) {
-    for (size_t j = 0; j < N; ++j) {
-      for (size_t k = 0; k < N; ++k) {
-
-        struct gpart* gp = &gparts[i * N * N + j * N + k];
-
-        gp->x[0] = i * boxSize / N + boxSize / (2 * N);
-        gp->x[1] = j * boxSize / N + boxSize / (2 * N);
-        gp->x[2] = k * boxSize / N + boxSize / (2 * N);
-
-        gp->mass = 1. / count;
-
-        gp->id_or_neg_offset = i * N * N + j * N + k;
-      }
-    }
-  }
-
-  /* Properties of the mesh */
-  const size_t meshmin[3] = {0, 0, 0};
-  const size_t meshmax[3] = {PMGRID - 1, PMGRID - 1, PMGRID - 1};
-
-  const size_t dimx = meshmax[0] - meshmin[0] + 2;
-  const size_t dimy = meshmax[1] - meshmin[1] + 2;
-  const size_t dimz = meshmax[2] - meshmin[2] + 2;
-
-  const double fac = PMGRID / boxSize;
-  const size_t PMGRID2 = 2 * (PMGRID / 2 + 1);
-
-  /* message("dimx=%zd dimy=%zd dimz=%zd", dimx, dimy, dimz); */
-
-  /* Allocate and empty the workspace mesh */
-  const size_t workspace_size = (dimx + 4) * (dimy + 4) * (dimz + 4);
-  double* workspace = fftw_malloc(workspace_size * sizeof(double));
-  bzero(workspace, workspace_size * sizeof(double));
-
-  /* Do CIC with the particles */
-  for (size_t pid = 0; pid < count; ++pid) {
-
-    const struct gpart* const gp = &gparts[pid];
-
-    const size_t slab_x =
-        (fac * gp->x[0] >= PMGRID) ? PMGRID - 1 : fac * gp->x[0];
-    const size_t slab_y =
-        (fac * gp->x[1] >= PMGRID) ? PMGRID - 1 : fac * gp->x[1];
-    const size_t slab_z =
-        (fac * gp->x[2] >= PMGRID) ? PMGRID - 1 : fac * gp->x[2];
-
-    const double dx = fac * gp->x[0] - (double)slab_x;
-    const double dy = fac * gp->x[1] - (double)slab_y;
-    const double dz = fac * gp->x[2] - (double)slab_z;
-
-    const size_t slab_xx = slab_x + 1;
-    const size_t slab_yy = slab_y + 1;
-    const size_t slab_zz = slab_z + 1;
-
-    workspace[(slab_x * dimy + slab_y) * dimz + slab_z] +=
-        gp->mass * (1.0 - dx) * (1.0 - dy) * (1.0 - dz);
-    workspace[(slab_x * dimy + slab_yy) * dimz + slab_z] +=
-        gp->mass * (1.0 - dx) * dy * (1.0 - dz);
-    workspace[(slab_x * dimy + slab_y) * dimz + slab_zz] +=
-        gp->mass * (1.0 - dx) * (1.0 - dy) * dz;
-    workspace[(slab_x * dimy + slab_yy) * dimz + slab_zz] +=
-        gp->mass * (1.0 - dx) * dy * dz;
-    workspace[(slab_xx * dimy + slab_y) * dimz + slab_z] +=
-        gp->mass * (dx) * (1.0 - dy) * (1.0 - dz);
-    workspace[(slab_xx * dimy + slab_yy) * dimz + slab_z] +=
-        gp->mass * (dx)*dy * (1.0 - dz);
-    workspace[(slab_xx * dimy + slab_y) * dimz + slab_zz] +=
-        gp->mass * (dx) * (1.0 - dy) * dz;
-    workspace[(slab_xx * dimy + slab_yy) * dimz + slab_zz] +=
-        gp->mass * (dx)*dy * dz;
-  }
-
-  /* for(size_t i = 0 ; i < dimx*dimy*dimz; ++i) */
-  /*   message("workspace[%zd] = %f", i, workspace[i]); */
-
-  /* Prepare the force grid */
-  const size_t fft_size = workspace_size;
-  double* forcegrid = fftw_malloc(fft_size * sizeof(double));
-  bzero(forcegrid, fft_size * sizeof(double));
-
-  const size_t sendmin = 0, recvmin = 0;
-  const size_t sendmax = PMGRID, recvmax = PMGRID;
-
-  memcpy(forcegrid, workspace + (sendmin - meshmin[0]) * dimy * dimz,
-         (sendmax - sendmin + 1) * dimy * dimz * sizeof(double));
-
-  /* for (size_t i = 0; i < fft_size; ++i) */
-  /*   if (forcegrid[i] != workspace[i]) error("wrong"); */
-
-  /* Prepare the density grid */
-  double* rhogrid = fftw_malloc(fft_size * sizeof(double));
-  bzero(rhogrid, fft_size * sizeof(double));
-
-  /* Now get the density */
-  for (size_t slab_x = recvmin; slab_x <= recvmax; slab_x++) {
-
-    const size_t slab_xx = slab_x % PMGRID;
-
-    for (size_t slab_y = recvmin; slab_y <= recvmax; slab_y++) {
-
-      const size_t slab_yy = slab_y % PMGRID;
-
-      for (size_t slab_z = recvmin; slab_z <= recvmax; slab_z++) {
-
-        const size_t slab_zz = slab_z % PMGRID;
-
-        rhogrid[PMGRID * PMGRID2 * slab_xx + PMGRID2 * slab_yy + slab_zz] +=
-            forcegrid[((slab_x - recvmin) * dimy + (slab_y - recvmin)) * dimz +
-                      (slab_z - recvmin)];
-      }
-    }
-  }
-
-  /* for (size_t i = 0; i < 640; i++) { */
-  /*   if (rhogrid[i] != workspace[i]) { */
-  /*     message("rhogrid[%zd]= %f workspace[%zd]= %f forcegrid[%zd]= %f", i, */
-  /*             rhogrid[i], i, workspace[i], i, forcegrid[i]); */
-  /*   } */
-  /* } */
-
-  /* FFT of the density field */
-  fftw_complex* fftgrid = fftw_malloc(fft_size * sizeof(fftw_complex));
-  fftw_plan plan_forward = fftw_plan_dft_r2c_3d(PMGRID, PMGRID, PMGRID, rhogrid,
-                                                fftgrid, FFTW_ESTIMATE);
-  fftw_execute(plan_forward);
-
-  for (size_t i = 0; i < 640; i++) {
-    message("workspace[%zd]= %f", i, fftgrid[i][0]);
+  /* Choke on FP-exceptions */
+  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+
+  /* Make one particle */
+  int nr_gparts = 1;
+  struct gpart *gparts = NULL;
+  posix_memalign((void **)&gparts, 64, nr_gparts * sizeof(struct gpart));
+  bzero(gparts, nr_gparts * sizeof(struct gpart));
+
+  gparts[0].x[0] = 0.3;
+  gparts[0].x[1] = 0.8;
+  gparts[0].x[2] = 0.2;
+  gparts[0].mass = 1.f;
+
+  /* Read the parameter file */
+  struct swift_params *params = malloc(sizeof(struct swift_params));
+  parser_read_file("fft_params.yml", params);
+
+  /* Initialise the gravity properties */
+  struct gravity_props gravity_properties;
+  gravity_props_init(&gravity_properties, params);
+
+  /* Build the infrastructure */
+  struct space space;
+  double dim[3] = {1., 1., 1.};
+  space_init(&space, params, dim, NULL, gparts, NULL, 0, nr_gparts, 0, 1, 1, 1,
+             0, 0);
+
+  struct engine engine;
+  engine.s = &space;
+  space.e = &engine;
+  engine.time = 0.1f;
+  engine.ti_current = 8;
+  engine.max_active_bin = num_time_bins;
+  engine.gravity_properties = &gravity_properties;
+  engine.nr_threads = 1;
+
+  struct runner runner;
+  runner.e = &engine;
+
+  /* Initialize the threadpool. */
+  threadpool_init(&engine.threadpool, engine.nr_threads);
+
+  space_rebuild(&space, 0);
+
+  /* Run the FFT task */
+  runner_do_grav_fft(&runner, 1);
+
+  /* Now check that we got the right answer */
+  int nr_cells = space.nr_cells;
+  double *r = malloc(nr_cells * sizeof(double));
+  double *pot = malloc(nr_cells * sizeof(double));
+  double *pot_exact = malloc(nr_cells * sizeof(double));
+
+  FILE *file = fopen("potential.dat", "w");
+  for (int i = 0; i < nr_cells; ++i) {
+    pot[i] = space.multipoles_top[i].pot.F_000;
+    double dx = NEAREST(space.multipoles_top[i].CoM[0] - gparts[0].x[0]);
+    double dy = NEAREST(space.multipoles_top[i].CoM[1] - gparts[0].x[1]);
+    double dz = NEAREST(space.multipoles_top[i].CoM[2] - gparts[0].x[2]);
+    r[i] = sqrt(dx * dx + dy * dy + dz * dz);
+    if (r[i] > 0) pot_exact[i] = -1. / r[i];
+    fprintf(file, "%e %e %e\n", r[i], pot[i], pot_exact[i]);
   }
+  fclose(file);
 
-  /* Clean-up */
-  fftw_destroy_plan(plan_forward);
-  fftw_free(forcegrid);
-  fftw_free(rhogrid);
-  fftw_free(workspace);
+  /* Clean up */
+  free(r);
+  free(pot);
+  free(pot_exact);
+  free(params);
   free(gparts);
   return 0;
 }