From dae880d619ec75d13577e7aed576d4106153c2bc Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Thu, 8 Mar 2018 16:25:52 +0900
Subject: [PATCH] Added calculation of the potential to the gravity exact
calculation for accuracy checks.
---
examples/main.c | 10 ++++----
src/gravity.c | 37 +++++++++++++++++++-----------
src/gravity/Default/gravity_part.h | 2 ++
src/kernel_gravity.h | 29 ++++++++++++++++++++---
4 files changed, 56 insertions(+), 22 deletions(-)
diff --git a/examples/main.c b/examples/main.c
index 828bf871bb..81211526df 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -802,6 +802,11 @@ int main(int argc, char *argv[]) {
e.dt_max);
fflush(stdout);
}
+
+/* Initialise the table of Ewald corrections for the gravity checks */
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ if (periodic) gravity_exact_force_ewald_init(e.s->dim[0]);
+#endif
}
/* Time to say good-bye if this was not a serious run. */
@@ -817,11 +822,6 @@ int main(int argc, char *argv[]) {
return 0;
}
-/* Initialise the table of Ewald corrections for the gravity checks */
-#ifdef SWIFT_GRAVITY_FORCE_CHECKS
- if (periodic) gravity_exact_force_ewald_init(e.s->dim[0]);
-#endif
-
/* Init the runner history. */
#ifdef HIST
for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
diff --git a/src/gravity.c b/src/gravity.c
index 05f4f37244..1e8ddc8c42 100644
--- a/src/gravity.c
+++ b/src/gravity.c
@@ -383,6 +383,7 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
/* Be ready for the calculation */
double a_grav[3] = {0., 0., 0.};
+ double pot = 0.;
/* Interact it with all other particles in the space.*/
for (int j = 0; j < (int)s->nr_gparts; ++j) {
@@ -408,27 +409,31 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
const double r_inv = 1. / sqrt(r2);
const double r = r2 * r_inv;
const double mj = gpj->mass;
- double f;
+ double f, phi;
if (r >= hi) {
/* Get Newtonian gravity */
f = mj * r_inv * r_inv * r_inv;
+ phi = -mj * r_inv;
} else {
const double ui = r * hi_inv;
- double W;
+ double Wf, Wp;
- kernel_grav_eval_double(ui, &W);
+ kernel_grav_eval_force_double(ui, &Wf);
+ kernel_grav_eval_pot_double(ui, &Wp);
/* Get softened gravity */
- f = mj * hi_inv3 * W;
+ f = mj * hi_inv3 * Wf;
+ phi = mj * hi_inv3 * Wp;
}
a_grav[0] += f * dx;
a_grav[1] += f * dy;
a_grav[2] += f * dz;
+ pot += phi;
/* Apply Ewald correction for periodic BC */
if (periodic && r > 1e-5 * hi) {
@@ -446,6 +451,7 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
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;
+ gpi->potential_exact = pot * const_G;
counter++;
}
@@ -529,8 +535,9 @@ 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\n", "id", "pos[0]",
- "pos[1]", "pos[2]", "a_swift[0]", "a_swift[1]", "a_swift[2]");
+ 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");
/* Output particle SWIFT accelerations */
for (size_t i = 0; i < s->nr_gparts; ++i) {
@@ -541,9 +548,10 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
if (gpi->id_or_neg_offset % 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 \n",
+ fprintf(file_swift,
+ "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n",
gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
- gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2]);
+ gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2], gpi->potential);
}
}
@@ -569,9 +577,9 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
fprintf(file_exact, "# periodic= %d\n", s->periodic);
fprintf(file_exact, "# Git Branch: %s\n", git_branch());
fprintf(file_exact, "# Git Revision: %s\n", git_revision());
- fprintf(file_exact, "# %16s %16s %16s %16s %16s %16s %16s\n", "id",
+ fprintf(file_exact, "# %16s %16s %16s %16s %16s %16s %16s %16s\n", "id",
"pos[0]", "pos[1]", "pos[2]", "a_exact[0]", "a_exact[1]",
- "a_exact[2]");
+ "a_exact[2]", "potential");
/* Output particle exact accelerations */
for (size_t i = 0; i < s->nr_gparts; ++i) {
@@ -582,10 +590,11 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
gpart_is_starting(gpi, e)) {
- fprintf(
- file_exact, "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e \n",
- gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
- gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2]);
+ fprintf(file_exact,
+ "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n",
+ gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
+ gpi->a_grav_exact[0], gpi->a_grav_exact[1],
+ gpi->a_grav_exact[2], gpi->potential_exact);
}
}
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 33c444e24e..4c7a7b5d4a 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -71,6 +71,8 @@ struct gpart {
/* Brute-force particle acceleration. */
double a_grav_exact[3];
+ /* Brute-force particle potential. */
+ double potential_exact;
#endif
} SWIFT_STRUCT_ALIGN;
diff --git a/src/kernel_gravity.h b/src/kernel_gravity.h
index aae694b8b2..24d30136c8 100644
--- a/src/kernel_gravity.h
+++ b/src/kernel_gravity.h
@@ -44,7 +44,7 @@ __attribute__((always_inline)) INLINE static void kernel_grav_pot_eval(
*W = *W * u;
*W = *W * u + 7.f;
*W = *W * u;
- *W = *W * u - 3;
+ *W = *W * u - 3.f;
}
/**
@@ -69,14 +69,37 @@ __attribute__((always_inline)) INLINE static void kernel_grav_force_eval(
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
/**
- * @brief Computes the gravity softening function in double precision.
+ * @brief Computes the gravity softening function for potential in double
+ * precision.
+ *
+ * This functions assumes 0 < u < 1.
+ *
+ * @param u The ratio of the distance to the softening length $u = x/h$.
+ * @param W (return) The value of the kernel function $W(x,h)$.
+ */
+__attribute__((always_inline)) INLINE static void kernel_grav_eval_pot_double(
+ double u, double *const W) {
+
+ /* W(u) = 3u^7 - 15u^6 + 28u^5 - 21u^4 + 7u^2 - 3 */
+ *W = 3. * u - 15.;
+ *W = *W * u + 28.;
+ *W = *W * u - 21.;
+ *W = *W * u;
+ *W = *W * u + 7.;
+ *W = *W * u;
+ *W = *W * u - 3;
+}
+
+/**
+ * @brief Computes the gravity softening function for forces in double
+ * precision.
*
* This functions assumes 0 < u < 1.
*
* @param u The ratio of the distance to the softening length $u = x/h$.
* @param W (return) The value of the kernel function $W(x,h)$.
*/
-__attribute__((always_inline)) INLINE static void kernel_grav_eval_double(
+__attribute__((always_inline)) INLINE static void kernel_grav_eval_force_double(
double u, double *const W) {
/* W(u) = 21u^5 - 90u^4 + 140u^3 - 84u^2 + 14 */
--
GitLab