diff --git a/src/gravity.c b/src/gravity.c
index 1642d37f7a6e543c6e2f263253aebc213ad4b19f..06a1e941670f1342db9780d9212b34b9f339f36c 100644
--- a/src/gravity.c
+++ b/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);
}
}
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index eabc95fc7eb1e1d4039dc9cea6611f1d17451a37..9f0db3f3bde9aef7a847d22dbc36b35b192b9304 100644
--- a/src/gravity/Default/gravity.h
+++ b/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
}
/**
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 06b9c52e14a655aa97b6c36ff9dca1f1cdd6e9a0..d325baf4de0938d8df539d38ae10f8f3f3ec7d2b 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/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];
diff --git a/src/runner_doiact_fft.c b/src/runner_doiact_fft.c
index 9bb5f799de5433f67bb3898291c6b55df31f1249..97283542f43a139afbff6073dde7a77cbcf094e2 100644
--- a/src/runner_doiact_fft.c
+++ b/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 */