Merge branch 'gravity_multi_dt' into 'master'

Gravity multi dt

This fully implements self-gravity bar the periodic task over MPI. 

Could you please check that is does not break anything on the hydro side ? 

If you are interested in running with self-gravity I'd recommend using the updated EAGLE_12 case and run with `-D -G`. 

My next steps is to work on the formal requirement of using `-D` alongside the `-G`. I also need to finish working out the automated scripts to generate code for the higher-order multipoles but that is orthogonal to the integration of the code into the master.

See merge request !324

configure.ac

@@ -799,6 +799,15 @@ case "$with_potential" in
    ;;
 esac
 
+#  Gravity multipole order
+AC_ARG_WITH([multipole-order],
+   [AS_HELP_STRING([--with-multipole-order=<order>],
+      [order of the multipole and gravitational field expansion @<:@ default: 3@:>@]
+   )],
+   [with_multipole_order="$withval"],
+   [with_multipole_order="3"]
+)
+AC_DEFINE_UNQUOTED([SELF_GRAVITY_MULTIPOLE_ORDER], [$with_multipole_order], [Multipole order])
 
 
 # Check for git, needed for revision stamps.
@@ -846,6 +855,7 @@ AC_MSG_RESULT([
    Riemann solver     : $with_riemann
    Cooling function   : $with_cooling
    External potential : $with_potential
+   Multipole order    : $with_multipole_order
 
    Task debugging     : $enable_task_debugging
    Debugging checks   : $enable_debugging_checks

examples/EAGLE_12/eagle_12.yml

@@ -13,6 +13,9 @@ TimeIntegration:
   dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-4  # The maximal time-step size of the simulation (in internal units).
 
+Scheduler:
+  cell_split_size:     50
+  
 # Parameters governing the snapshots
 Snapshots:
   basename:            eagle # Common part of the name of output files
@@ -23,6 +26,13 @@ Snapshots:
 Statistics:
   delta_time:          1e-2 # Time between statistics output
 
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                   0.025    # Constant dimensionless multiplier for time integration. 
+  epsilon:               0.0001   # Softening length (in internal units).
+  a_smooth:              1000.
+  r_cut:                 4.
+  
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).

examples/EAGLE_25/eagle_25.yml

@@ -23,6 +23,13 @@ Snapshots:
 Statistics:
   delta_time:          1e-2 # Time between statistics output
 
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                   0.025    # Constant dimensionless multiplier for time integration. 
+  epsilon:               0.0001   # Softening length (in internal units).
+  a_smooth:              1000.
+  r_cut:                 4.
+  
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).

examples/UniformDMBox/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.985  ,
+'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-6
+max_error = 1e-1
+num_bins = 51
+
+# 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 = ['b', 'g', 'r', 'm']
+
+# Time-step to plot
+step = int(sys.argv[1])
+
+# Find the files for the different expansion orders
+order_list = glob.glob("gravity_checks_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][26])
+    
+# Start the plot
+plt.figure()
+
+# 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, :]
+    
+    # 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_per95 = np.percentile(norm_error,95)
+    per95_x = np.percentile(error_x,95)
+    per95_y = np.percentile(error_y,95)
+    per95_z = np.percentile(error_z,95)
+
+    plt.subplot(221)    
+    plt.semilogx(bins, norm_error_hist, 'k--', label="Gadget-2")
+    plt.plot([norm_median, norm_median], [2.7, 3], 'k-', lw=1)
+    plt.plot([norm_per95, norm_per95], [2.7, 3], 'k:', lw=1)
+    plt.subplot(222)
+    plt.semilogx(bins, error_x_hist, 'k--', label="Gadget-2")
+    plt.plot([median_x, median_x], [1.8, 2], 'k-', lw=1)
+    plt.plot([per95_x, per95_x], [1.8, 2], 'k:', lw=1)
+    plt.subplot(223)    
+    plt.semilogx(bins, error_y_hist, 'k--', label="Gadget-2")
+    plt.plot([median_y, median_y], [1.8, 2], 'k-', lw=1)
+    plt.plot([per95_y, per95_y], [1.8, 2], 'k:', lw=1)
+    plt.subplot(224)    
+    plt.semilogx(bins, error_z_hist, 'k--', label="Gadget-2")
+    plt.plot([median_z, median_z], [1.8, 2], 'k-', lw=1)
+    plt.plot([per95_z, per95_z], [1.8, 2], 'k:', lw=1)
+
+
+# Plot the different histograms
+for i in range(num_order-1, -1, -1):
+    data = np.loadtxt(order_list[i])
+    ids = data[:,0]
+    pos = data[:,1:4]
+    a_exact = data[:,4:7]
+    a_grav = data[:, 7:10]
+
+    # Sort stuff
+    sort_index = np.argsort(ids)
+    ids = ids[sort_index]
+    pos = pos[sort_index, :]
+    a_exact = a_exact[sort_index, :]
+    a_grav = a_grav[sort_index, :]
+
+    # Cross-checks
+    if not np.array_equal(ids, gadget2_ids):
+        print "Comparing different IDs !"
+
+    if not np.array_equal(pos, gadget2_pos):
+        print "Comparing different positions ! max difference:", np.max(pos - gadget2_pos)
+
+    if not np.array_equal(a_exact, gadget2_a_exact):
+        print "Comparing different exact accelerations ! max difference:", np.max(a_exact - gadget2_a_exact)
+        
+        
+    # Compute the error norm
+    diff = a_exact - a_grav
+
+    norm_diff = np.sqrt(diff[:,0]**2 + diff[:,1]**2 + diff[:,2]**2)
+    norm_a = np.sqrt(a_exact[:,0]**2 + a_exact[:,1]**2 + 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_per95 = np.percentile(norm_error,95)
+    per95_x = np.percentile(error_x,95)
+    per95_y = np.percentile(error_y,95)
+    per95_z = np.percentile(error_z,95)
+    
+    plt.subplot(221)
+    plt.semilogx(bins, norm_error_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+    plt.plot([norm_median, norm_median], [2.7, 3],'-', color=cols[i], lw=1)
+    plt.plot([norm_per95, norm_per95], [2.7, 3],':', color=cols[i], lw=1)
+    plt.subplot(222)    
+    plt.semilogx(bins, error_x_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+    plt.plot([median_x, median_x], [1.8, 2],'-', color=cols[i], lw=1)
+    plt.plot([per95_x, per95_x], [1.8, 2],':', color=cols[i], lw=1)
+    plt.subplot(223)    
+    plt.semilogx(bins, error_y_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+    plt.plot([median_y, median_y], [1.8, 2],'-', color=cols[i], lw=1)
+    plt.plot([per95_y, per95_y], [1.8, 2],':', color=cols[i], lw=1)
+    plt.subplot(224)    
+    plt.semilogx(bins, error_z_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+    plt.plot([median_z, median_z], [1.8, 2],'-', color=cols[i], lw=1)
+    plt.plot([per95_z, per95_z], [1.8, 2],':', color=cols[i], lw=1)
+
+    
+plt.subplot(221)
+plt.xlabel("$|\delta \overrightarrow{a}|/|\overrightarrow{a}_{exact}|$")
+plt.ylabel("Density")
+plt.xlim(min_error, 2*max_error)
+plt.ylim(0,3)
+plt.legend(loc="upper left")
+plt.subplot(222)    
+plt.xlabel("$\delta a_x/|\overrightarrow{a}_{exact}|$")
+plt.ylabel("Density")
+plt.xlim(min_error, 2*max_error)
+plt.ylim(0,2)
+#plt.legend(loc="center left")
+plt.subplot(223)    
+plt.xlabel("$\delta a_y/|\overrightarrow{a}_{exact}|$")
+plt.ylabel("Density")
+plt.xlim(min_error, 2*max_error)
+plt.ylim(0,2)
+#plt.legend(loc="center left")
+plt.subplot(224)    
+plt.xlabel("$\delta a_z/|\overrightarrow{a}_{exact}|$")
+plt.ylabel("Density")
+plt.xlim(min_error, 2*max_error)
+plt.ylim(0,2)
+#plt.legend(loc="center left")
+
+
+
+plt.savefig("gravity_checks_step%d.png"%step)

examples/main.c

@@ -323,14 +323,14 @@ int main(int argc, char *argv[]) {
 
   /* How large are the parts? */
   if (myrank == 0) {
-    message("sizeof(part)       is %4zi bytes.", sizeof(struct part));
-    message("sizeof(xpart)      is %4zi bytes.", sizeof(struct xpart));
-    message("sizeof(spart)      is %4zi bytes.", sizeof(struct spart));
-    message("sizeof(gpart)      is %4zi bytes.", sizeof(struct gpart));
-    message("sizeof(multipole)  is %4zi bytes.", sizeof(struct multipole));
-    message("sizeof(acc_tensor) is %4zi bytes.", sizeof(struct acc_tensor));
-    message("sizeof(task)       is %4zi bytes.", sizeof(struct task));
-    message("sizeof(cell)       is %4zi bytes.", sizeof(struct cell));
+    message("sizeof(part)        is %4zi bytes.", sizeof(struct part));
+    message("sizeof(xpart)       is %4zi bytes.", sizeof(struct xpart));
+    message("sizeof(spart)       is %4zi bytes.", sizeof(struct spart));
+    message("sizeof(gpart)       is %4zi bytes.", sizeof(struct gpart));
+    message("sizeof(multipole)   is %4zi bytes.", sizeof(struct multipole));
+    message("sizeof(grav_tensor) is %4zi bytes.", sizeof(struct grav_tensor));
+    message("sizeof(task)        is %4zi bytes.", sizeof(struct task));
+    message("sizeof(cell)        is %4zi bytes.", sizeof(struct cell));
   }
 
   /* Read the parameter file */

examples/plot_tasks.py

@@ -57,14 +57,15 @@ pl.rcParams.update(PLOT_PARAMS)
 #  Tasks and subtypes. Indexed as in tasks.h.
 TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
              "init", "ghost", "extra_ghost", "drift", "kick1", "kick2",
-             "timestep", "send", "recv", "grav_gather_m", "grav_fft",
-             "grav_mm", "grav_up", "cooling", "sourceterms", "count"]
+             "timestep", "send", "recv", "grav_top_level", "grav_long_range",
+             "grav_mm", "grav_down", "cooling", "sourceterms", "count"]
 SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
-            "tend", "xv", "rho", "gpart", "count"]
+            "tend", "xv", "rho", "gpart", "multipole", "spart", "count"]
 
 #  Task/subtypes of interest.
-FULLTYPES = ["self/force", "self/density", "sub_self/force",
-             "sub_self/density", "pair/force", "pair/density", "sub_pair/force",
+FULLTYPES = ["self/force", "self/density", "self/grav", "sub_self/force",
+             "sub_self/density", "pair/force", "pair/density", "pair/grav",
+             "sub_pair/force",
              "sub_pair/density", "recv/xv", "send/xv", "recv/rho", "send/rho",
              "recv/tend", "send/tend"]
 

src/Makefile.am

@@ -46,7 +46,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
     sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
     dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
-    hydro_space.h
+    vector_power.h hydro_space.h
 
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \

src/cell.c

@@ -1114,7 +1114,7 @@ void cell_check_multipole(struct cell *c, void *data) {
 
 #ifdef SWIFT_DEBUG_CHECKS
   struct gravity_tensors ma;
-  const double tolerance = 1e-5; /* Relative */
+  const double tolerance = 1e-3; /* Relative */
 
   /* First recurse */
   if (c->split)
@@ -1127,8 +1127,7 @@ void cell_check_multipole(struct cell *c, void *data) {
     gravity_P2M(&ma, c->gparts, c->gcount);
 
     /* Now  compare the multipole expansion */
-    if (!gravity_multipole_equal(&ma.m_pole, &c->multipole->m_pole,
-                                 tolerance)) {
+    if (!gravity_multipole_equal(&ma, c->multipole, tolerance)) {
       message("Multipoles are not equal at depth=%d!", c->depth);
       message("Correct answer:");
       gravity_multipole_print(&ma.m_pole);
@@ -1487,7 +1486,7 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
   } else if (ti_current > ti_old_multipole) {
 
     /* Drift the multipole */
-    gravity_multipole_drift(c->multipole, dt);
+    gravity_drift(c->multipole, dt);
   }
 
   /* Update the time of the last drift */
@@ -1504,7 +1503,21 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
  * @param e The #engine (to get ti_current).
  */
 void cell_drift_multipole(struct cell *c, const struct engine *e) {
-  error("To be implemented");
+
+  const double timeBase = e->timeBase;
+  const integertime_t ti_old_multipole = c->ti_old_multipole;
+  const integertime_t ti_current = e->ti_current;
+
+  /* Drift from the last time the cell was drifted to the current time */
+  const double dt = (ti_current - ti_old_multipole) * timeBase;
+
+  /* Check that we are actually going to move forward. */
+  if (ti_current < ti_old_multipole) error("Attempt to drift to the past");
+
+  if (ti_current > ti_old_multipole) gravity_drift(c->multipole, dt);
+
+  /* Update the time of the last drift */
+  c->ti_old_multipole = ti_current;
 }
 
 /**

src/const.h

@@ -39,9 +39,6 @@
 /* Thermal energy per unit mass used as a constant for the isothermal EoS */
 #define const_isothermal_internal_energy 20.2615290634f
 
-/* Self gravity stuff. */
-#define const_gravity_multipole_order 1
-
 /* Type of gradients to use (GIZMO_SPH only) */
 /* If no option is chosen, no gradients are used (first order scheme) */
 //#define GRADIENTS_SPH

src/engine.c

@@ -3029,6 +3029,13 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
 
   if (e->nodeID == 0) message("Computing initial gas densities.");
 
+  /* Initialise the softening lengths */
+  if (e->policy & engine_policy_self_gravity) {
+
+    for (size_t i = 0; i < s->nr_gparts; ++i)
+      gravity_init_softening(&s->gparts[i], e->gravity_properties);
+  }
+
   /* Construct all cells and tasks to start everything */
   engine_rebuild(e);
 
@@ -3061,16 +3068,17 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
   }
 
 #ifdef SWIFT_DEBUG_CHECKS
-  /* Let's store the total mass in the system for future checks */
-  e->s->total_mass = 0.;
-  for (size_t i = 0; i < s->nr_gparts; ++i)
-    e->s->total_mass += s->gparts[i].mass;
-#ifdef WITH_MPI
-  if (MPI_Allreduce(MPI_IN_PLACE, &e->s->total_mass, 1, MPI_DOUBLE, MPI_SUM,
-                    MPI_COMM_WORLD) != MPI_SUCCESS)
-    error("Failed to all-reduce total mass in the system.");
-#endif
-  message("Total mass in the system: %e", e->s->total_mass);
+  /* Check that we have the correct total mass in the top-level multipoles */
+  size_t num_gpart_mpole = 0;
+  if (e->policy & engine_policy_self_gravity) {
+    for (int i = 0; i < e->s->nr_cells; ++i)
+      num_gpart_mpole += e->s->cells_top[i].multipole->m_pole.num_gpart;
+    if (num_gpart_mpole != e->s->nr_gparts)
+      error(
+          "Multipoles don't contain the total number of gpart s->nr_gpart=%zd, "
+          "m_poles=%zd",
+          e->s->nr_gparts, num_gpart_mpole);
+  }
 #endif
 
   /* Now time to get ready for the first time-step */
@@ -3085,9 +3093,19 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
   /* Print the number of active tasks ? */
   if (e->verbose) engine_print_task_counts(e);
 
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  /* Run the brute-force gravity calculation for some gparts */
+  gravity_exact_force_compute(e->s, e);
+#endif
+
   /* Run the 0th time-step */
   engine_launch(e, e->nr_threads);
 
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  /* Check the accuracy of the gravity calculation */
+  gravity_exact_force_check(e->s, e, 1e-1);
+#endif
+
   /* Recover the (integer) end of the next time-step */
   engine_collect_timestep(e);
 
@@ -3160,6 +3178,17 @@ void engine_step(struct engine *e) {
   /* Print the number of active tasks ? */
   if (e->verbose) engine_print_task_counts(e);
 
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Check that we have the correct total mass in the top-level multipoles */
+  size_t num_gpart_mpole = 0;
+  if (e->policy & engine_policy_self_gravity) {
+    for (int i = 0; i < e->s->nr_cells; ++i)
+      num_gpart_mpole += e->s->cells_top[i].multipole->m_pole.num_gpart;
+    if (num_gpart_mpole != e->s->nr_gparts)
+      error("Multipoles don't contain the total number of gpart");
+  }
+#endif
+
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   /* Run the brute-force gravity calculation for some gparts */
   gravity_exact_force_compute(e->s, e);

src/gravity.c

@@ -20,6 +20,9 @@
 /* Config parameters. */
 #include "../config.h"
 
+/* Some standard headers. */
+#include <stdio.h>
+
 /* This object's header. */
 #include "gravity.h"
 
@@ -53,9 +56,7 @@ void gravity_exact_force_compute(struct space *s, const struct engine *e) {
         gpart_is_active(gpi, e)) {
 
       /* Be ready for the calculation */
-      gpi->a_grav[0] = 0.f;
-      gpi->a_grav[1] = 0.f;
-      gpi->a_grav[2] = 0.f;
+      double a_grav[3] = {0., 0., 0.};
 
       /* Interact it with all other particles in the space.*/
       for (size_t j = 0; j < s->nr_gparts; ++j) {
@@ -66,36 +67,55 @@ void gravity_exact_force_compute(struct space *s, const struct engine *e) {
         struct gpart *gpj = &s->gparts[j];
 
         /* Compute the pairwise distance. */
-        float dx[3] = {gpi->x[0] - gpj->x[0],   // x
-                       gpi->x[1] - gpj->x[1],   // y
-                       gpi->x[2] - gpj->x[2]};  // z
-        const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+        const double dx[3] = {gpi->x[0] - gpj->x[0],   // x
+                              gpi->x[1] - gpj->x[1],   // y
+                              gpi->x[2] - gpj->x[2]};  // z
+        const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
 
-        runner_iact_grav_pp_nonsym(0.f, r2, dx, gpi, gpj);
-      }
+        const double r = sqrtf(r2);
+        const double ir = 1.f / r;
+        const double mj = gpj->mass;
+        const double hi = gpi->epsilon;
+        double f;
+        const double f_lr = 1.;
 
-      /* Finish the calculation */
-      gravity_end_force(gpi, const_G);
+        if (r >= hi) {
 
-      /* Store the exact answer */
-      gpi->a_grav_exact[0] = gpi->a_grav[0];
-      gpi->a_grav_exact[1] = gpi->a_grav[1];
-      gpi->a_grav_exact[2] = gpi->a_grav[2];
+          /* Get Newtonian gravity */
+          f = mj * ir * ir * ir * f_lr;
+
+        } else {
+
+          const double hi_inv = 1.f / hi;
+          const double hi_inv3 = hi_inv * hi_inv * hi_inv;
+          const double ui = r * hi_inv;
+          float W;
+
+          kernel_grav_eval(ui, &W);
 
-      /* Restore everything */
-      gpi->a_grav[0] = 0.f;
-      gpi->a_grav[1] = 0.f;
-      gpi->a_grav[2] = 0.f;
+          /* Get softened gravity */
+          f = mj * hi_inv3 * W * f_lr;
+        }
+
+        const double fdx[3] = {f * dx[0], f * dx[1], f * dx[2]};
+
+        a_grav[0] -= fdx[0];
+        a_grav[1] -= fdx[1];
+        a_grav[2] -= fdx[2];
+      }
+
+      /* Store the exact answer */
+      gpi->a_grav_exact[0] = a_grav[0] * const_G;
+      gpi->a_grav_exact[1] = a_grav[1] * const_G;
+      gpi->a_grav_exact[2] = a_grav[2] * const_G;
 
       counter++;
     }
   }
 
-  message("Computed exact gravity for %d gparts.", counter);
+  message("Computed exact gravity for %d gparts (took %.3f %s). ", counter,
+          clocks_from_ticks(getticks() - tic), clocks_getunit());
 
-  if (e->verbose)
-    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
-            clocks_getunit());
 #else
   error("Gravity checking function called without the corresponding flag.");
 #endif
@@ -118,17 +138,24 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
 
-  const double const_G = e->physical_constants->const_newton_G;
-
   int counter = 0;
 
   /* Some accumulators */
-  float err_rel[3];
-  float err_rel_max[3] = {0.f, 0.f, 0.f};
-  float err_rel_min[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
-  float err_rel_mean[3] = {0.f, 0.f, 0.f};
-  float err_rel_mean2[3] = {0.f, 0.f, 0.f};
-  float err_rel_std[3] = {0.f, 0.f, 0.f};
+  float err_rel_max = 0.f;
+  float err_rel_min = FLT_MAX;
+  float err_rel_mean = 0.f;
+  float err_rel_mean2 = 0.f;
+  float err_rel_std = 0.f;
+