Added direct calculation of the mesh forces when checking accuracy. Use...

Added direct calculation of the mesh forces when checking accuracy. Use higher-order CIC interpolation when translating the mesh back to the field tensors.

src/gravity.c

@@ -588,9 +588,11 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
   fprintf(file_swift, "# theta= %16.8e\n", e->gravity_properties->theta_crit);
   fprintf(file_swift, "# Git Branch: %s\n", git_branch());
   fprintf(file_swift, "# Git Revision: %s\n", git_revision());
-  fprintf(file_swift, "# %16s %16s %16s %16s %16s %16s %16s %16s\n", "id",
-          "pos[0]", "pos[1]", "pos[2]", "a_swift[0]", "a_swift[1]",
-          "a_swift[2]", "potential");
+  fprintf(file_swift,
+          "# %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s\n",
+          "id", "pos[0]", "pos[1]", "pos[2]", "a_swift[0]", "a_swift[1]",
+          "a_swift[2]", "potential", "a_PM[0]", "a_PM[1]", "a_PM[2]",
+          "potentialPM");
 
   /* Output particle SWIFT accelerations  */
   for (size_t i = 0; i < s->nr_gparts; ++i) {
@@ -611,9 +613,11 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
     if (id % SWIFT_GRAVITY_FORCE_CHECKS == 0 && gpart_is_starting(gpi, e)) {
 
       fprintf(file_swift,
-              "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n", id,
-              gpi->x[0], gpi->x[1], gpi->x[2], gpi->a_grav[0], gpi->a_grav[1],
-              gpi->a_grav[2], gpi->potential);
+              "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e "
+              "%16.8e %16.8e %16.8e\n",
+              id, gpi->x[0], gpi->x[1], gpi->x[2], gpi->a_grav[0],
+              gpi->a_grav[1], gpi->a_grav[2], gpi->potential, gpi->a_grav_PM[0],
+              gpi->a_grav_PM[1], gpi->a_grav_PM[2], gpi->potential_PM);
     }
   }
 

src/gravity/Default/gravity.h

@@ -130,6 +130,13 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
   gp->a_grav[2] = 0.f;
   gp->potential = 0.f;
 
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  gp->potential_PM = 0.f;
+  gp->a_grav_PM[0] = 0.f;
+  gp->a_grav_PM[1] = 0.f;
+  gp->a_grav_PM[2] = 0.f;
+#endif
+
 #ifdef SWIFT_DEBUG_CHECKS
   gp->num_interacted = 0;
 #endif
@@ -151,6 +158,13 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
   gp->a_grav[1] *= const_G;
   gp->a_grav[2] *= const_G;
   gp->potential *= const_G;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  gp->potential_PM *= const_G;
+  gp->a_grav_PM[0] *= const_G;
+  gp->a_grav_PM[1] *= const_G;
+  gp->a_grav_PM[2] *= const_G;
+#endif
 }
 
 /**

src/gravity/Default/gravity_part.h

@@ -65,6 +65,12 @@ struct gpart {
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
 
+  /*! Acceleration taken from the mesh only */
+  float a_grav_PM[3];
+
+  /*! Potential taken from the mesh only */
+  float potential_PM;
+
   /* Brute-force particle acceleration. */
   double a_grav_exact[3];
 

src/runner_doiact_fft.c

@@ -28,9 +28,11 @@
 #include "runner_doiact_fft.h"
 
 /* Local includes. */
+#include "debug.h"
 #include "engine.h"
 #include "error.h"
 #include "kernel_long_gravity.h"
+#include "part.h"
 #include "runner.h"
 #include "space.h"
 #include "timers.h"
@@ -40,6 +42,9 @@
 /**
  * @brief Returns 1D index of a 3D NxNxN array using row-major style.
  *
+ * Wraps around in the corresponding dimension if any of the 3 indices is >= N
+ * or < 0.
+ *
  * @param i Index along x.
  * @param j Index along y.
  * @param k Index along z.
@@ -47,7 +52,38 @@
  */
 __attribute__((always_inline)) INLINE static int row_major_id(int i, int j,
                                                               int k, int N) {
-  return ((i % N) * N * N + (j % N) * N + (k % N));
+  return (((i + N) % N) * N * N + ((j + N) % N) * N + ((k + N) % N));
+}
+
+/**
+ * @brief Interpolate values from a the mesh using CIC.
+ *
+ * @param mesh The mesh to read from.
+ * @param i The index of the cell along x
+ * @param j The index of the cell along y
+ * @param k The index of the cell along z
+ * @param tx First CIC coefficient along x
+ * @param ty First CIC coefficient along y
+ * @param tz First CIC coefficient along z
+ * @param dx Second CIC coefficient along x
+ * @param dy Second CIC coefficient along y
+ * @param dz Second CIC coefficient along z
+ */
+__attribute__((always_inline)) INLINE static double CIC_get(
+    double mesh[6][6][6], int i, int j, int k, double tx, double ty, double tz,
+    double dx, double dy, double dz) {
+
+  double temp;
+  temp = mesh[i + 0][j + 0][k + 0] * tx * ty * tz;
+  temp += mesh[i + 0][j + 0][k + 1] * tx * ty * dz;
+  temp += mesh[i + 0][j + 1][k + 0] * tx * dy * tz;
+  temp += mesh[i + 0][j + 1][k + 1] * tx * dy * dz;
+  temp += mesh[i + 1][j + 0][k + 0] * dx * ty * tz;
+  temp += mesh[i + 1][j + 0][k + 1] * dx * ty * dz;
+  temp += mesh[i + 1][j + 1][k + 0] * dx * dy * tz;
+  temp += mesh[i + 1][j + 1][k + 1] * dx * dy * dz;
+
+  return temp;
 }
 
 /**
@@ -89,17 +125,18 @@ __attribute__((always_inline)) INLINE static void multipole_to_mesh_CIC(
   if (k < 0 || k >= N) error("Invalid multipole position in z");
 #endif
 
+  const double mass = m->m_pole.M_000;
+
   /* CIC ! */
-  rho[row_major_id(i + 0, j + 0, k + 0, N)] += m->m_pole.M_000 * tx * ty * tz;
-  rho[row_major_id(i + 0, j + 0, k + 1, N)] += m->m_pole.M_000 * tx * ty * dz;
-  rho[row_major_id(i + 0, j + 1, k + 0, N)] += m->m_pole.M_000 * tx * dy * tz;
-  rho[row_major_id(i + 0, j + 1, k + 1, N)] += m->m_pole.M_000 * tx * dy * dz;
-  rho[row_major_id(i + 1, j + 0, k + 0, N)] += m->m_pole.M_000 * dx * ty * tz;
-  rho[row_major_id(i + 1, j + 0, k + 1, N)] += m->m_pole.M_000 * dx * ty * dz;
-  rho[row_major_id(i + 1, j + 1, k + 0, N)] += m->m_pole.M_000 * dx * dy * tz;
-  rho[row_major_id(i + 1, j + 1, k + 1, N)] += m->m_pole.M_000 * dx * dy * dz;
+  rho[row_major_id(i + 0, j + 0, k + 0, N)] += mass * tx * ty * tz;
+  rho[row_major_id(i + 0, j + 0, k + 1, N)] += mass * tx * ty * dz;
+  rho[row_major_id(i + 0, j + 1, k + 0, N)] += mass * tx * dy * tz;
+  rho[row_major_id(i + 0, j + 1, k + 1, N)] += mass * tx * dy * dz;
+  rho[row_major_id(i + 1, j + 0, k + 0, N)] += mass * dx * ty * tz;
+  rho[row_major_id(i + 1, j + 0, k + 1, N)] += mass * dx * ty * dz;
+  rho[row_major_id(i + 1, j + 1, k + 0, N)] += mass * dx * dy * tz;
+  rho[row_major_id(i + 1, j + 1, k + 1, N)] += mass * dx * dy * dz;
 }
-
 /**
  * @brief Computes the potential on a multipole from a given mesh using the CIC
  * method.
@@ -140,19 +177,263 @@ __attribute__((always_inline)) INLINE static void mesh_to_multipole_CIC(
   if (k < 0 || k >= N) error("Invalid multipole position in z");
 #endif
 
-  /* CIC ! */
-  m->pot.F_000 += pot[row_major_id(i + 0, j + 0, k + 0, N)] * tx * ty * tz;
-  m->pot.F_000 += pot[row_major_id(i + 0, j + 0, k + 1, N)] * tx * ty * dz;
-  m->pot.F_000 += pot[row_major_id(i + 0, j + 1, k + 0, N)] * tx * dy * tz;
-  m->pot.F_000 += pot[row_major_id(i + 0, j + 1, k + 1, N)] * tx * dy * dz;
-  m->pot.F_000 += pot[row_major_id(i + 1, j + 0, k + 0, N)] * dx * ty * tz;
-  m->pot.F_000 += pot[row_major_id(i + 1, j + 0, k + 1, N)] * dx * ty * dz;
-  m->pot.F_000 += pot[row_major_id(i + 1, j + 1, k + 0, N)] * dx * dy * tz;
-  m->pot.F_000 += pot[row_major_id(i + 1, j + 1, k + 1, N)] * dx * dy * dz;
+  /* First, copy the necessary part of the mesh for stencil operations */
+  /* This includes box-wrapping in all 3 dimensions. */
+  double phi[6][6][6];
+  for (int iii = -2; iii <= 3; ++iii) {
+    for (int jjj = -2; jjj <= 3; ++jjj) {
+      for (int kkk = -2; kkk <= 3; ++kkk) {
+        phi[iii + 2][jjj + 2][kkk + 2] =
+            pot[row_major_id(i + iii, j + jjj, k + kkk, N)];
+      }
+    }
+  }
+
+  /* Indices of (i,j,k) in the local copy of the mesh */
+  const int ii = 2, jj = 2, kk = 2;
+
+  /* Some local accumulators */
+  double F_000 = 0.;
+  double F_100 = 0., F_010 = 0., F_001 = 0.;
+  double F_200 = 0., F_020 = 0., F_002 = 0.;
+  double F_110 = 0., F_101 = 0., F_011 = 0.;
+  double F_300 = 0., F_030 = 0., F_003 = 0.;
+
+  /* Simple CIC for the potential itself */
+  F_000 -= CIC_get(phi, ii, jj, kk, tx, ty, tz, dx, dy, dz);
+
+  /* ---- */
+
+  /* 5-point stencil along each axis for the 1st derivatives */
+  F_100 -= (1. / 12.) * CIC_get(phi, ii + 2, jj, kk, tx, ty, tz, dx, dy, dz);
+  F_100 += (2. / 3.) * CIC_get(phi, ii + 1, jj, kk, tx, ty, tz, dx, dy, dz);
+  F_100 -= (2. / 3.) * CIC_get(phi, ii - 1, jj, kk, tx, ty, tz, dx, dy, dz);
+  F_100 += (1. / 12.) * CIC_get(phi, ii - 2, jj, kk, tx, ty, tz, dx, dy, dz);
+
+  F_010 -= (1. / 12.) * CIC_get(phi, ii, jj + 2, kk, tx, ty, tz, dx, dy, dz);
+  F_010 += (2. / 3.) * CIC_get(phi, ii, jj + 1, kk, tx, ty, tz, dx, dy, dz);
+  F_010 -= (2. / 3.) * CIC_get(phi, ii, jj - 1, kk, tx, ty, tz, dx, dy, dz);
+  F_010 += (1. / 12.) * CIC_get(phi, ii, jj - 2, kk, tx, ty, tz, dx, dy, dz);
+
+  F_001 -= (1. / 12.) * CIC_get(phi, ii, jj, kk + 2, tx, ty, tz, dx, dy, dz);
+  F_001 += (2. / 3.) * CIC_get(phi, ii, jj, kk + 1, tx, ty, tz, dx, dy, dz);
+  F_001 -= (2. / 3.) * CIC_get(phi, ii, jj, kk - 1, tx, ty, tz, dx, dy, dz);
+  F_001 += (1. / 12.) * CIC_get(phi, ii, jj, kk - 2, tx, ty, tz, dx, dy, dz);
+
+  /* ---- */
+
+  /* 5-point stencil along each axis for the 2nd derivatives (diagonal) */
+  F_200 -= (1. / 12.) * CIC_get(phi, ii + 2, jj, kk, tx, ty, tz, dx, dy, dz);
+  F_200 += (4. / 3.) * CIC_get(phi, ii + 1, jj, kk, tx, ty, tz, dx, dy, dz);
+  F_200 -= (5. / 2.) * CIC_get(phi, ii + 0, jj, kk, tx, ty, tz, dx, dy, dz);
+  F_200 += (4. / 3.) * CIC_get(phi, ii - 1, jj, kk, tx, ty, tz, dx, dy, dz);
+  F_200 -= (1. / 12.) * CIC_get(phi, ii - 2, jj, kk, tx, ty, tz, dx, dy, dz);
+
+  F_020 -= (1. / 12.) * CIC_get(phi, ii, jj + 2, kk, tx, ty, tz, dx, dy, dz);
+  F_020 += (4. / 3.) * CIC_get(phi, ii, jj + 1, kk, tx, ty, tz, dx, dy, dz);
+  F_020 -= (5. / 2.) * CIC_get(phi, ii, jj + 0, kk, tx, ty, tz, dx, dy, dz);
+  F_020 += (4. / 3.) * CIC_get(phi, ii, jj - 1, kk, tx, ty, tz, dx, dy, dz);
+  F_020 -= (1. / 12.) * CIC_get(phi, ii, jj - 2, kk, tx, ty, tz, dx, dy, dz);
+
+  F_002 -= (1. / 12.) * CIC_get(phi, ii, jj, kk + 2, tx, ty, tz, dx, dy, dz);
+  F_002 += (4. / 3.) * CIC_get(phi, ii, jj, kk + 1, tx, ty, tz, dx, dy, dz);
+  F_002 -= (5. / 2.) * CIC_get(phi, ii, jj, kk + 0, tx, ty, tz, dx, dy, dz);
+  F_002 += (4. / 3.) * CIC_get(phi, ii, jj, kk - 1, tx, ty, tz, dx, dy, dz);
+  F_002 -= (1. / 12.) * CIC_get(phi, ii, jj, kk - 2, tx, ty, tz, dx, dy, dz);
+
+  /* Regular stencil for the 2nd derivatives (non-diagonal) */
+  F_110 += (1. / 4.) * CIC_get(phi, ii + 1, jj + 1, kk, tx, ty, tz, dx, dy, dz);
+  F_110 -= (1. / 4.) * CIC_get(phi, ii + 1, jj - 1, kk, tx, ty, tz, dx, dy, dz);
+  F_110 -= (1. / 4.) * CIC_get(phi, ii - 1, jj + 1, kk, tx, ty, tz, dx, dy, dz);
+  F_110 += (1. / 4.) * CIC_get(phi, ii - 1, jj - 1, kk, tx, ty, tz, dx, dy, dz);
+
+  F_101 += (1. / 4.) * CIC_get(phi, ii + 1, jj, kk + 1, tx, ty, tz, dx, dy, dz);
+  F_101 -= (1. / 4.) * CIC_get(phi, ii + 1, jj, kk - 1, tx, ty, tz, dx, dy, dz);
+  F_101 -= (1. / 4.) * CIC_get(phi, ii - 1, jj, kk + 1, tx, ty, tz, dx, dy, dz);
+  F_101 += (1. / 4.) * CIC_get(phi, ii - 1, jj, kk - 1, tx, ty, tz, dx, dy, dz);
+
+  F_011 += (1. / 4.) * CIC_get(phi, ii, jj + 1, kk + 1, tx, ty, tz, dx, dy, dz);
+  F_011 -= (1. / 4.) * CIC_get(phi, ii, jj + 1, kk - 1, tx, ty, tz, dx, dy, dz);
+  F_011 -= (1. / 4.) * CIC_get(phi, ii, jj - 1, kk + 1, tx, ty, tz, dx, dy, dz);
+  F_011 += (1. / 4.) * CIC_get(phi, ii, jj - 1, kk - 1, tx, ty, tz, dx, dy, dz);
+
+  /* ---- */
+
+  /* 5-point stencil along each axis for the 3rd derivatives (diagonal) */
+  F_300 -= (1. / 2.) * CIC_get(phi, ii + 2, jj, kk, tx, ty, tz, dx, dy, dz);
+  F_300 += 1. * CIC_get(phi, ii + 1, jj, kk, tx, ty, tz, dx, dy, dz);
+  F_300 -= 1. * CIC_get(phi, ii - 1, jj, kk, tx, ty, tz, dx, dy, dz);
+  F_300 += (1. / 2.) * CIC_get(phi, ii - 2, jj, kk, tx, ty, tz, dx, dy, dz);
+
+  F_030 -= (1. / 2.) * CIC_get(phi, ii, jj + 2, kk, tx, ty, tz, dx, dy, dz);
+  F_030 += (2. / 3.) * CIC_get(phi, ii, jj + 1, kk, tx, ty, tz, dx, dy, dz);
+  F_030 -= (2. / 3.) * CIC_get(phi, ii, jj - 1, kk, tx, ty, tz, dx, dy, dz);
+  F_030 += (1. / 2.) * CIC_get(phi, ii, jj - 2, kk, tx, ty, tz, dx, dy, dz);
+
+  F_003 -= (1. / 2.) * CIC_get(phi, ii, jj, kk + 2, tx, ty, tz, dx, dy, dz);
+  F_003 += 1. * CIC_get(phi, ii, jj, kk + 1, tx, ty, tz, dx, dy, dz);
+  F_003 -= 1. * CIC_get(phi, ii, jj, kk - 1, tx, ty, tz, dx, dy, dz);
+  F_003 += (1. / 2.) * CIC_get(phi, ii, jj, kk - 2, tx, ty, tz, dx, dy, dz);
+
+  /* Store things back */
+  m->pot.F_000 += F_000;
+  m->pot.F_100 -= F_100 * fac;
+  m->pot.F_010 -= F_010 * fac;
+  m->pot.F_001 -= F_001 * fac;
+  m->pot.F_200 += F_200 * fac * fac;
+  m->pot.F_020 += F_020 * fac * fac;
+  m->pot.F_002 += F_002 * fac * fac;
+  m->pot.F_110 -= F_110 * fac * fac;
+  m->pot.F_011 -= F_011 * fac * fac;
+  m->pot.F_101 -= F_101 * fac * fac;
+  m->pot.F_300 += F_300 * fac * fac * fac;
+  m->pot.F_030 += F_030 * fac * fac * fac;
+  m->pot.F_003 += F_003 * fac * fac * fac;
+
+  m->pot.interacted = 1;
 }
 
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+
+/**
+ * @brief Computes the potential on a gpart from a given mesh using the CIC
+ * method.
+ *
+ * Debugging routine.
+ *
+ * @param gp The #gpart.
+ * @param pot The potential mesh.
+ * @param N the size of the mesh along one axis.
+ * @param fac width of a mesh cell.
+ * @param dim The dimensions of the simulation box.
+ */
+void mesh_to_gparts_CIC(struct gpart* gp, const double* pot, int N, double fac,
+                        const double dim[3]) {
+
+  /* Box wrap the multipole's position */
+  const double pos_x = box_wrap(gp->x[0], 0., dim[0]);
+  const double pos_y = box_wrap(gp->x[1], 0., dim[1]);
+  const double pos_z = box_wrap(gp->x[2], 0., dim[2]);
+
+  int i = (int)(fac * pos_x);
+  if (i >= N) i = N - 1;
+  const double dx = fac * pos_x - i;
+  const double tx = 1. - dx;
+
+  int j = (int)(fac * pos_y);
+  if (j >= N) j = N - 1;
+  const double dy = fac * pos_y - j;
+  const double ty = 1. - dy;
+
+  int k = (int)(fac * pos_z);
+  if (k >= N) k = N - 1;
+  const double dz = fac * pos_z - k;
+  const double tz = 1. - dz;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (i < 0 || i >= N) error("Invalid multipole position in x");
+  if (j < 0 || j >= N) error("Invalid multipole position in y");
+  if (k < 0 || k >= N) error("Invalid multipole position in z");
 #endif
 
+  if (gp->a_grav_PM[0] != 0. || gp->potential_PM != 0.)
+    error("Particle with non-initalised stuff");
+
+  /* First, copy the necessary part of the mesh for stencil operations */
+  /* This includes box-wrapping in all 3 dimensions. */
+  double phi[6][6][6];
+  for (int iii = -2; iii <= 3; ++iii) {
+    for (int jjj = -2; jjj <= 3; ++jjj) {
+      for (int kkk = -2; kkk <= 3; ++kkk) {
+        phi[iii + 2][jjj + 2][kkk + 2] =
+            pot[row_major_id(i + iii, j + jjj, k + kkk, N)];
+      }
+    }
+  }
+
+  /* Some local accumulators */
+  double p = 0.;
+  double a[3] = {0.};
+
+  /* Indices of (i,j,k) in the local copy of the mesh */
+  const int ii = 2, jj = 2, kk = 2;
+
+  /* Simple CIC for the potential itself */
+  p += CIC_get(phi, ii, jj, kk, tx, ty, tz, dx, dy, dz);
+
+  /* ---- */
+
+  /* 5-point stencil along each axis for the accelerations */
+  a[0] += (1. / 12.) * CIC_get(phi, ii + 2, jj, kk, tx, ty, tz, dx, dy, dz);
+  a[0] -= (2. / 3.) * CIC_get(phi, ii + 1, jj, kk, tx, ty, tz, dx, dy, dz);
+  a[0] += (2. / 3.) * CIC_get(phi, ii - 1, jj, kk, tx, ty, tz, dx, dy, dz);
+  a[0] -= (1. / 12.) * CIC_get(phi, ii - 2, jj, kk, tx, ty, tz, dx, dy, dz);
+
+  a[1] += (1. / 12.) * CIC_get(phi, ii, jj + 2, kk, tx, ty, tz, dx, dy, dz);
+  a[1] -= (2. / 3.) * CIC_get(phi, ii, jj + 1, kk, tx, ty, tz, dx, dy, dz);
+  a[1] += (2. / 3.) * CIC_get(phi, ii, jj - 1, kk, tx, ty, tz, dx, dy, dz);
+  a[1] -= (1. / 12.) * CIC_get(phi, ii, jj - 2, kk, tx, ty, tz, dx, dy, dz);
+
+  a[2] += (1. / 12.) * CIC_get(phi, ii, jj, kk + 2, tx, ty, tz, dx, dy, dz);
+  a[2] -= (2. / 3.) * CIC_get(phi, ii, jj, kk + 1, tx, ty, tz, dx, dy, dz);
+  a[2] += (2. / 3.) * CIC_get(phi, ii, jj, kk - 1, tx, ty, tz, dx, dy, dz);
+  a[2] -= (1. / 12.) * CIC_get(phi, ii, jj, kk - 2, tx, ty, tz, dx, dy, dz);
+
+  /* ---- */
+
+  /* Store things back */
+  gp->potential_PM = p;
+  gp->a_grav_PM[0] = fac * a[0];
+  gp->a_grav_PM[1] = fac * a[1];
+  gp->a_grav_PM[2] = fac * a[2];
+}
+#endif
+
+/**
+ * @brief Dump a real array of size NxNxN to stdout.
+ *
+ * Debugging routine.
+ *
+ * @param array The array to dump.
+ * @param N The side-length of the array to dump
+ */
+void print_array(double* array, int N) {
+
+  for (int k = N - 1; k >= 0; --k) {
+    printf("--- z = %d ---------\n", k);
+    for (int j = N - 1; j >= 0; --j) {
+      for (int i = 0; i < N; ++i) {
+        printf("%f ", array[i * N * N + j * N + k]);
+      }
+      printf("\n");
+    }
+  }
+}
+
+/**
+ * @brief Dump a complex array of size NxNxN to stdout.
+ *
+ * Debugging routine.
+ *
+ * @param array The array to dump.
+ * @param N The side-length of the array to dump
+ */
+void print_carray(fftw_complex* array, int N) {
+
+  for (int k = N - 1; k >= 0; --k) {
+    printf("--- z = %d ---------\n", k);
+    for (int j = N - 1; j >= 0; --j) {
+      for (int i = 0; i < N; ++i) {
+        printf("(%f %f) ", array[i * N * N + j * N + k][0],
+               array[i * N * N + j * N + k][1]);
+      }
+      printf("\n");
+    }
+  }
+}
+
+#endif /* HAVE_FFTW */
+
 /**
  * @brief Computes the potential on the top multipoles using a Fourier transform
  *
@@ -214,6 +495,8 @@ void runner_do_grav_fft(struct runner* r, int timer) {
   for (int i = 0; i < nr_cells; ++i)
     multipole_to_mesh_CIC(&multipoles[i], rho, N, cell_fac, dim);
 
+  /* print_array(rho, N); */
+
   /* Fourier transform to go to magic-land */
   fftw_execute(forward_plan);
 
@@ -286,11 +569,22 @@ void runner_do_grav_fft(struct runner* r, int timer) {
   fftw_execute(inverse_plan);
 
   /* rho now contains the potential */
+  /* Let's create an alias to avoid confusion */
   /* This array is now again NxNxN real numbers */
+  double* potential = rho;
+
+  /* message("\n\n\n POTENTIAL"); */
+  /* print_array(potential, N); */
 
-  /* Get the potential from the mesh using CIC */
+  /* Get the potential from the mesh to the gravity tensors using CIC */
   for (int i = 0; i < nr_cells; ++i)
-    mesh_to_multipole_CIC(&multipoles[i], rho, N, cell_fac, dim);
+    mesh_to_multipole_CIC(&multipoles[i], potential, N, cell_fac, dim);
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  /* Get the potential from the mesh to the gparts using CIC */
+  for (size_t i = 0; i < s->nr_gparts; ++i)
+    mesh_to_gparts_CIC(&s->gparts[i], potential, N, cell_fac, dim);
+#endif
 
   /* Clean-up the mess */
   fftw_destroy_plan(forward_plan);
@@ -305,32 +599,3 @@ void runner_do_grav_fft(struct runner* r, int timer) {
   error("No FFTW library found. Cannot compute periodic long-range forces.");
 #endif
 }
-
-#ifdef HAVE_FFTW
-void print_array(double* array, int N) {
-
-  for (int k = N - 1; k >= 0; --k) {
-    printf("--- z = %d ---------\n", k);
-    for (int j = N - 1; j >= 0; --j) {
-      for (int i = 0; i < N; ++i) {
-        printf("%f ", array[i * N * N + j * N + k]);
-      }
-      printf("\n");
-    }
-  }
-}
-
-void print_carray(fftw_complex* array, int N) {
-
-  for (int k = N - 1; k >= 0; --k) {
-    printf("--- z = %d ---------\n", k);
-    for (int j = N - 1; j >= 0; --j) {
-      for (int i = 0; i < N; ++i) {
-        printf("(%f %f) ", array[i * N * N + j * N + k][0],
-               array[i * N * N + j * N + k][1]);
-      }
-      printf("\n");
-    }
-  }
-}
-#endif /* HAVE_FFTW */