Merge branch 'gravity_multi_dt' into 'master'

Gravity multi dt

Some improvements to the gravity code:

 - Corrected typo in 4th order vector powers.
 - Use a user-defined opening angle for the distance checks.
 - Increase accuracy of the exact gravity calculation.
 - Only do exact gravity calculation if a matching result file is not present in the same directory.
 - Reinstate the dumping of individual task timers.
 - Added script to plot the accuracy of the calculation.
 - Crash if the M2L kernel is called with a length smaller than the softening length.

See merge request !331

.gitignore

@@ -32,6 +32,7 @@ examples/*/*/*.hdf5
 examples/*/*/*.png
 examples/*/*/*.txt
 examples/*/*/used_parameters.yml
+examples/*/gravity_checks_*.dat
 
 tests/testPair
 tests/brute_force_standard.dat

README

@@ -31,6 +31,7 @@ Valid options are:
   -s          Run with hydrodynamics.
   -S          Run with stars.
   -t    {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified.
+  -T          Print timers every time-step.
   -v     [12] Increase the level of verbosity
   	      1: MPI-rank 0 writes
 	      2: All MPI-ranks write

examples/EAGLE_12/eagle_12.yml

@@ -28,7 +28,8 @@ Statistics:
 
 # Parameters for the self-gravity scheme
 Gravity:
-  eta:                   0.025    # Constant dimensionless multiplier for time integration. 
+  eta:                   0.025    # Constant dimensionless multiplier for time integration.
+  theta:                 0.7      # Opening angle (Multipole acceptance criterion)
   epsilon:               0.0001   # Softening length (in internal units).
   a_smooth:              1000.
   r_cut:                 4.

examples/UniformDMBox/uniformBox.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.    # The end time of the simulation (in internal units).
+  dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-3  # The maximal time-step size of the simulation (in internal units).
+
+Scheduler:
+  max_top_level_cells: 8
+  cell_split_size:     50
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            uniformDMBox # Common part of the name of output files
+  time_first:          0.           # Time of the first output (in internal units)
+  delta_time:          0.01         # Time difference between consecutive outputs (in internal units)
+
+# 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).
+  
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./uniformDMBox_100.hdf5     # The file to read

examples/main.c

@@ -90,6 +90,7 @@ void print_help_message() {
   printf("  %2s %8s %s\n", "-t", "{int}",
          "The number of threads to use on each MPI rank. Defaults to 1 if not "
          "specified.");
+  printf("  %2s %8s %s\n", "-T", "", "Print timers every time-step.");
   printf("  %2s %8s %s\n", "-v", "[12]", "Increase the level of verbosity.");
   printf("  %2s %8s %s\n", "", "", "1: MPI-rank 0 writes ");
   printf("  %2s %8s %s\n", "", "", "2: All MPI-ranks write");
@@ -165,12 +166,13 @@ int main(int argc, char *argv[]) {
   int with_drift_all = 0;
   int verbose = 0;
   int nr_threads = 1;
+  int with_verbose_timers = 0;
   char paramFileName[200] = "";
   unsigned long long cpufreq = 0;
 
   /* Parse the parameters */
   int c;
-  while ((c = getopt(argc, argv, "acCdDef:FgGhn:sSt:v:y:")) != -1) switch (c) {
+  while ((c = getopt(argc, argv, "acCdDef:FgGhn:sSt:Tv:y:")) != -1) switch (c) {
       case 'a':
         with_aff = 1;
         break;
@@ -229,6 +231,9 @@ int main(int argc, char *argv[]) {
           return 1;
         }
         break;
+      case 'T':
+        with_verbose_timers = 1;
+        break;
       case 'v':
         if (sscanf(optarg, "%d", &verbose) != 1) {
           if (myrank == 0) printf("Error parsing verbosity level (-v).\n");
@@ -596,11 +601,18 @@ int main(int argc, char *argv[]) {
   engine_dump_snapshot(&e);
 
   /* Legend */
-  if (myrank == 0)
+  if (myrank == 0) {
     printf("# %6s %14s %14s %10s %10s %10s %16s [%s]\n", "Step", "Time",
            "Time-step", "Updates", "g-Updates", "s-Updates", "Wall-clock time",
            clocks_getunit());
 
+    if (with_verbose_timers) {
+      printf("timers: ");
+      for (int k = 0; k < timer_count; k++) printf("%s\t", timers_names[k]);
+      printf("\n");
+    }
+  }
+
   /* Main simulation loop */
   for (int j = 0; !engine_is_done(&e) && e.step - 1 != nsteps; j++) {
 
@@ -610,6 +622,15 @@ int main(int argc, char *argv[]) {
     /* Take a step. */
     engine_step(&e);
 
+    /* Print the timers. */
+    if (with_verbose_timers) {
+      printf("timers: ");
+      for (int k = 0; k < timer_count; k++)
+        printf("%.3f\t", clocks_from_ticks(timers[k]));
+      printf("\n");
+      timers_reset(0xFFFFFFFFllu);
+    }
+
 #ifdef SWIFT_DEBUG_TASKS
     /* Dump the task data using the given frequency. */
     if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {

examples/parameter_example.yml

@@ -52,6 +52,7 @@ SPH:
 # 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.1      # Softening length (in internal units).
   a_smooth:              1.25     # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
   r_cut:                 4.5      # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).

examples/plot_gravity_checks.py

+#!/usr/bin/env python
+
+import sys
+import glob
+import re
+import numpy as np
+import matplotlib.pyplot as plt
+
+params = {'axes.labelsize': 14,
+'axes.titlesize': 18,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 14,
+'ytick.labelsize': 14,
+'text.usetex': True,
+'figure.figsize': (10, 10),
+'figure.subplot.left'    : 0.06,
+'figure.subplot.right'   : 0.99  ,
+'figure.subplot.bottom'  : 0.06  ,
+'figure.subplot.top'     : 0.99  ,
+'figure.subplot.wspace'  : 0.14  ,
+'figure.subplot.hspace'  : 0.14  ,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+plt.rcParams.update(params)
+plt.rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+min_error = 1e-7
+max_error = 3e-1
+num_bins = 64
+
+# Construct the bins
+bin_edges = np.linspace(np.log10(min_error), np.log10(max_error), num_bins + 1)
+bin_size = (np.log10(max_error) - np.log10(min_error)) / num_bins
+bins = 0.5*(bin_edges[1:] + bin_edges[:-1])
+bin_edges = 10**bin_edges
+bins = 10**bins
+
+# Colours
+cols = ['#332288', '#88CCEE', '#117733', '#DDCC77', '#CC6677']
+
+# Time-step to plot
+step = int(sys.argv[1])
+
+# Find the files for the different expansion orders
+order_list = glob.glob("gravity_checks_swift_step%d_order*.dat"%step)
+num_order = len(order_list)
+
+# Get the multipole orders
+order = np.zeros(num_order)
+for i in range(num_order):
+    order[i] = int(order_list[i][32])
+order = sorted(order)
+order_list = sorted(order_list)
+
+# Read the exact accelerations first
+data = np.loadtxt('gravity_checks_exact_step%d.dat'%step)
+exact_ids = data[:,0]
+exact_pos = data[:,1:4]
+exact_a = data[:,4:7]
+# Sort stuff
+sort_index = np.argsort(exact_ids)
+exact_ids = exact_ids[sort_index]
+exact_pos = exact_pos[sort_index, :]
+exact_a = exact_a[sort_index, :]        
+exact_a_norm = np.sqrt(exact_a[:,0]**2 + exact_a[:,1]**2 + exact_a[:,2]**2)
+    
+# Start the plot
+plt.figure()
+
+count = 0
+
+# Get the Gadget-2 data if existing
+gadget2_file_list = glob.glob("forcetest_gadget2.txt")
+if len(gadget2_file_list) != 0:
+
+    gadget2_data = np.loadtxt(gadget2_file_list[0])
+    gadget2_ids = gadget2_data[:,0]
+    gadget2_pos = gadget2_data[:,1:4]
+    gadget2_a_exact = gadget2_data[:,4:7]
+    gadget2_a_grav = gadget2_data[:, 7:10]
+
+    # Sort stuff
+    sort_index = np.argsort(gadget2_ids)
+    gadget2_ids = gadget2_ids[sort_index]
+    gadget2_pos = gadget2_pos[sort_index, :]
+    gadget2_a_exact = gadget2_a_exact[sort_index, :]
+    gadget2_a_grav = gadget2_a_grav[sort_index, :]
+
+    # Cross-checks
+    if not np.array_equal(exact_ids, gadget2_ids):
+        print "Comparing different IDs !"
+
+    if np.max(np.abs(exact_pos - gadget2_pos)/np.abs(gadget2_pos)) > 1e-6:
+        print "Comparing different positions ! max difference:"
+        index = np.argmax(exact_pos[:,0]**2 + exact_pos[:,1]**2 + exact_pos[:,2]**2 - gadget2_pos[:,0]**2 - gadget2_pos[:,1]**2 - gadget2_pos[:,2]**2)
+        print "Gadget2 (id=%d):"%gadget2_ids[index], gadget2_pos[index,:], "exact (id=%d):"%exact_ids[index], exact_pos[index,:], "\n"
+
+    if np.max(np.abs(exact_a - gadget2_a_exact) / np.abs(gadget2_a_exact)) > 2e-6:
+        print "Comparing different exact accelerations ! max difference:"
+        index = np.argmax(exact_a[:,0]**2 + exact_a[:,1]**2 + exact_a[:,2]**2 - gadget2_a_exact[:,0]**2 - gadget2_a_exact[:,1]**2 - gadget2_a_exact[:,2]**2)
+        print "a_exact --- Gadget2:", gadget2_a_exact[index,:], "exact:", exact_a[index,:]
+        print "pos ---     Gadget2: (id=%d):"%gadget2_ids[index], gadget2_pos[index,:], "exact (id=%d):"%ids[index], pos[index,:],"\n"
+
+    
+    # Compute the error norm
+    diff = gadget2_a_exact - gadget2_a_grav
+
+    norm_diff = np.sqrt(diff[:,0]**2 + diff[:,1]**2 + diff[:,2]**2)
+    norm_a = np.sqrt(gadget2_a_exact[:,0]**2 + gadget2_a_exact[:,1]**2 + gadget2_a_exact[:,2]**2)
+
+    norm_error = norm_diff / norm_a
+    error_x = abs(diff[:,0]) / norm_a
+    error_y = abs(diff[:,1]) / norm_a
+    error_z = abs(diff[:,2]) / norm_a
+    
+    # Bin the error
+    norm_error_hist,_ = np.histogram(norm_error, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+    error_x_hist,_ = np.histogram(error_x, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+    error_y_hist,_ = np.histogram(error_y, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+    error_z_hist,_ = np.histogram(error_z, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+    
+    norm_median = np.median(norm_error)
+    median_x = np.median(error_x)
+    median_y = np.median(error_y)
+    median_z = np.median(error_z)
+
+    norm_per99 = np.percentile(norm_error,99)
+    per99_x = np.percentile(error_x,99)
+    per99_y = np.percentile(error_y,99)
+    per99_z = np.percentile(error_z,99)
+
+    norm_max = np.max(norm_error)
+    max_x = np.max(error_x)
+    max_y = np.max(error_y)
+    max_z = np.max(error_z)
+
+    print "Gadget-2 ---- "
+    print "Norm: median= %f 99%%= %f max= %f"%(norm_median, norm_per99, norm_max)
+    print "X   : median= %f 99%%= %f max= %f"%(median_x, per99_x, max_x)
+    print "Y   : median= %f 99%%= %f max= %f"%(median_y, per99_y, max_y)
+    print "Z   : median= %f 99%%= %f max= %f"%(median_z, per99_z, max_z)
+    print ""
+
+    plt.subplot(221)    
+    plt.text(min_error * 1.5, 1.55, "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(norm_median, norm_per99), ha="left", va="top", alpha=0.8)
+    plt.semilogx(bins, norm_error_hist, 'k--', label="Gadget-2", alpha=0.8)
+    plt.subplot(222)
+    plt.semilogx(bins, error_x_hist, 'k--', label="Gadget-2", alpha=0.8)
+    plt.text(min_error * 1.5, 1.55, "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_x, per99_x), ha="left", va="top", alpha=0.8)
+    plt.subplot(223)    
+    plt.semilogx(bins, error_y_hist, 'k--', label="Gadget-2", alpha=0.8)
+    plt.text(min_error * 1.5, 1.55, "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_y, per99_y), ha="left", va="top", alpha=0.8)
+    plt.subplot(224)    
+    plt.semilogx(bins, error_z_hist, 'k--', label="Gadget-2", alpha=0.8)
+    plt.text(min_error * 1.5, 1.55, "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_z, per99_z), ha="left", va="top", alpha=0.8)
+    
+    count += 1
+
+
+# Plot the different histograms
+for i in range(num_order):
+    data = np.loadtxt(order_list[i])
+    ids = data[:,0]
+    pos = data[:,1:4]
+    a_grav = data[:, 4:7]
+
+    # Sort stuff
+    sort_index = np.argsort(ids)
+    ids = ids[sort_index]
+    pos = pos[sort_index, :]
+    a_grav = a_grav[sort_index, :]        
+
+    # Cross-checks
+    if not np.array_equal(exact_ids, ids):
+        print "Comparing different IDs !"
+
+    if np.max(np.abs(exact_pos - pos)/np.abs(pos)) > 1e-6:
+        print "Comparing different positions ! max difference:"
+        index = np.argmax(exact_pos[:,0]**2 + exact_pos[:,1]**2 + exact_pos[:,2]**2 - pos[:,0]**2 - pos[:,1]**2 - pos[:,2]**2)
+        print "SWIFT (id=%d):"%ids[index], pos[index,:], "exact (id=%d):"%exact_ids[index], exact_pos[index,:], "\n"
+
+    
+    # Compute the error norm
+    diff = exact_a - a_grav
+
+    norm_diff = np.sqrt(diff[:,0]**2 + diff[:,1]**2 + diff[:,2]**2)
+
+    norm_error = norm_diff / exact_a_norm
+    error_x = abs(diff[:,0]) / exact_a_norm
+    error_y = abs(diff[:,1]) / exact_a_norm
+    error_z = abs(diff[:,2]) / exact_a_norm
+    
+    # Bin the error
+    norm_error_hist,_ = np.histogram(norm_error, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+    error_x_hist,_ = np.histogram(error_x, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+    error_y_hist,_ = np.histogram(error_y, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+    error_z_hist,_ = np.histogram(error_z, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+
+    norm_median = np.median(norm_error)
+    median_x = np.median(error_x)
+    median_y = np.median(error_y)
+    median_z = np.median(error_z)
+
+    norm_per99 = np.percentile(norm_error,99)
+    per99_x = np.percentile(error_x,99)
+    per99_y = np.percentile(error_y,99)
+    per99_z = np.percentile(error_z,99)
+
+    norm_max = np.max(norm_error)
+    max_x = np.max(error_x)
+    max_y = np.max(error_y)
+    max_z = np.max(error_z)
+
+    print "Order %d ---- "%order[i]
+    print "Norm: median= %f 99%%= %f max= %f"%(norm_median, norm_per99, norm_max)
+    print "X   : median= %f 99%%= %f max= %f"%(median_x, per99_x, max_x)
+    print "Y   : median= %f 99%%= %f max= %f"%(median_y, per99_y, max_y)
+    print "Z   : median= %f 99%%= %f max= %f"%(median_z, per99_z, max_z)
+    print ""
+    
+    plt.subplot(221)
+    plt.semilogx(bins, norm_error_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+    plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(norm_median, norm_per99), ha="left", va="top", color=cols[i])
+    plt.subplot(222)    
+    plt.semilogx(bins, error_x_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+    plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_x, per99_x), ha="left", va="top", color=cols[i])
+    plt.subplot(223)    
+    plt.semilogx(bins, error_y_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+    plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_y, per99_y), ha="left", va="top", color=cols[i])
+    plt.subplot(224)    
+    plt.semilogx(bins, error_z_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+    plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_z, per99_z), ha="left", va="top", color=cols[i])
+
+    count += 1
+
+plt.subplot(221)
+plt.xlabel("$|\delta \overrightarrow{a}|/|\overrightarrow{a}_{exact}|$")
+#plt.ylabel("Density")
+plt.xlim(min_error, max_error)
+plt.ylim(0,2.5)
+plt.legend(loc="upper left")
+plt.subplot(222)    
+plt.xlabel("$\delta a_x/|\overrightarrow{a}_{exact}|$")
+#plt.ylabel("Density")
+plt.xlim(min_error, max_error)
+plt.ylim(0,1.75)
+#plt.legend(loc="center left")
+plt.subplot(223)    
+plt.xlabel("$\delta a_y/|\overrightarrow{a}_{exact}|$")
+#plt.ylabel("Density")
+plt.xlim(min_error, max_error)
+plt.ylim(0,1.75)
+#plt.legend(loc="center left")
+plt.subplot(224)    
+plt.xlabel("$\delta a_z/|\overrightarrow{a}_{exact}|$")
+#plt.ylabel("Density")
+plt.xlim(min_error, max_error)
+plt.ylim(0,1.75)
+#plt.legend(loc="center left")
+
+
+
+plt.savefig("gravity_checks_step%d.png"%step, dpi=200)
+plt.savefig("gravity_checks_step%d.pdf"%step, dpi=200)

src/cell.c

@@ -1165,6 +1165,16 @@ void cell_check_multipole(struct cell *c, void *data) {
       gravity_multipole_print(&c->multipole->m_pole);
       error("Aborting");
     }
+
+    /* Check that the upper limit of r_max is good enough */
+    if (!(c->multipole->r_max >= ma.r_max)) {
+      error("Upper-limit r_max=%e too small. Should be >=%e.",
+            c->multipole->r_max, ma.r_max);
+    } else if (c->multipole->r_max * c->multipole->r_max >
+               3. * c->width[0] * c->width[0]) {
+      error("r_max=%e larger than cell diagonal %e.", c->multipole->r_max,
+            sqrt(3. * c->width[0] * c->width[0]));
+    }
   }
 #else
   error("Calling debugging code without debugging flag activated.");

src/engine.c

@@ -1647,9 +1647,9 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
  * @brief Constructs the top-level tasks for the short-range gravity
  * interactions.
  *
- * All top-cells get a self task.
- * All neighbouring pairs get a pair task.
- * All non-neighbouring pairs within a range of 6 cells get a M-M task.
+ * - All top-cells get a self task.
+ * - All pairs within range according to the multipole acceptance
+ *   criterion get a pair task.
  *
  * @param e The #engine.
  */
@@ -1658,6 +1658,7 @@ void engine_make_self_gravity_tasks(struct engine *e) {
   struct space *s = e->s;
   struct scheduler *sched = &e->sched;
   const int nodeID = e->nodeID;
+  const double theta_crit_inv = e->gravity_properties->theta_crit_inv;
   struct cell *cells = s->cells_top;
   const int nr_cells = s->nr_cells;
 
@@ -1668,13 +1669,14 @@ void engine_make_self_gravity_tasks(struct engine *e) {
     /* Skip cells without gravity particles */
     if (ci->gcount == 0) continue;
 
-    /* Is that neighbour local ? */
+    /* Is that cell local ? */
     if (ci->nodeID != nodeID) continue;
 
     /* If the cells is local build a self-interaction */
     scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci, NULL,
                       0);
 
+    /* Loop over every other cell */
     for (int cjd = cid + 1; cjd < nr_cells; ++cjd) {
 
       struct cell *cj = &cells[cjd];
@@ -1683,9 +1685,11 @@ void engine_make_self_gravity_tasks(struct engine *e) {
       if (cj->gcount == 0) continue;
 
       /* Is that neighbour local ? */
-      if (cj->nodeID != nodeID) continue;
+      if (cj->nodeID != nodeID) continue;  // MATTHIEU
 
-      if (cell_are_neighbours(ci, cj))
+      /* Are the cells to close for a MM interaction ? */
+      if (!gravity_multipole_accept(ci->multipole, cj->multipole,
+                                    theta_crit_inv))
         scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, ci,
                           cj, 1);
     }

src/gravity.c

@@ -22,6 +22,7 @@
 
 /* Some standard headers. */
 #include <stdio.h>
+#include <unistd.h>
 
 /* This object's header. */
 #include "gravity.h"
@@ -29,6 +30,57 @@
 /* Local headers. */
 #include "active.h"
 #include "error.h"
+#include "version.h"
+
+/**
+ * @brief Checks whether the file containing the exact accelerations for
+ * the current choice of parameters already exists.
+ *
+ * @param e The #engine.
+ */
+int gravity_exact_force_file_exits(const struct engine *e) {
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+
+  /* File name */
+  char file_name[100];
+  sprintf(file_name, "gravity_checks_exact_step%d.dat", e->step);
+
+  /* Does the file exist ? */
+  if (access(file_name, R_OK | W_OK) == 0) {
+
+    /* Let's check whether the header matches the parameters of this run */
+    FILE *file = fopen(file_name, "r");
+    if (!file) error("Problem reading gravity_check file");
+
+    char line[100];
+    char dummy1[10], dummy2[10];
+    double epsilon, newton_G;
+    int N;
+    /* Reads file header */
+    if (fgets(line, 100, file) != line) error("Problem reading title");
+    if (fgets(line, 100, file) != line) error("Problem reading G");
+    sscanf(line, "%s %s %le", dummy1, dummy2, &newton_G);
+    if (fgets(line, 100, file) != line) error("Problem reading N");
+    sscanf(line, "%s %s %d", dummy1, dummy2, &N);
+    if (fgets(line, 100, file) != line) error("Problem reading epsilon");
+    sscanf(line, "%s %s %le", dummy1, dummy2, &epsilon);
+    fclose(file);
+
+    /* Check whether it matches the current parameters */
+    if (N == SWIFT_GRAVITY_FORCE_CHECKS &&
+        (fabs(epsilon - e->gravity_properties->epsilon) / epsilon < 1e-5) &&
+        (fabs(newton_G - e->physical_constants->const_newton_G) / newton_G <
+         1e-5)) {
+      return 1;
+    }
+  }
+  return 0;
+#else
+  error("Gravity checking function called without the corresponding flag.");
+  return 0;
+#endif
+}
 
 /**
  * @brief Run a brute-force gravity calculation for a subset of particles.
@@ -47,6 +99,13 @@ void gravity_exact_force_compute(struct space *s, const struct engine *e) {
   const double const_G = e->physical_constants->const_newton_G;
   int counter = 0;
 
+  /* Let's start by checking whether we already computed these forces */
+  if (gravity_exact_force_file_exits(e)) {
+    message("Exact accelerations already computed. Skipping calculation.");
+    return;
+  }
+
+  /* No matching file present ? Do it then */
   for (size_t i = 0; i < s->nr_gparts; ++i) {
 
     struct gpart *gpi = &s->gparts[i];
@@ -72,8 +131,8 @@ void gravity_exact_force_compute(struct space *s, const struct engine *e) {
                               gpi->x[2] - gpj->x[2]};  // z
         const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
 
-        const double r = sqrtf(r2);
-        const double ir = 1.f / r;
+        const double r = sqrt(r2);
+        const double ir = 1. / r;
         const double mj = gpj->mass;
         const double hi = gpi->epsilon;
         double f;
@@ -86,15 +145,17 @@ void gravity_exact_force_compute(struct space *s, const struct engine *e) {
 
         } else {
 
-          const double hi_inv = 1.f / hi;
+          const double hi_inv = 1. / hi;
           const double hi_inv3 = hi_inv * hi_inv * hi_inv;
           const double ui = r * hi_inv;
-          float W;
+          double W;
 
-          kernel_grav_eval(ui, &W);
+          kernel_grav_eval_double(ui, &W);
 
           /* Get softened gravity */
           f = mj * hi_inv3 * W * f_lr;
+
+          // printf("r=%e hi=%e W=%e fac=%e\n", r, hi, W, f);
         }
 
         const double fdx[3] = {f * dx[0], f * dx[1], f * dx[2]};
@@ -138,25 +199,24 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
 
-  int counter = 0;
-
-  /* Some accumulators */
-  float err_rel_max = 0.f;