Added check for gravity accuracy for a user-specified fraction of all particles

examples/main.c

@@ -305,6 +305,15 @@ int main(int argc, char *argv[]) {
     message("WARNING: Debugging checks activated. Code will be slower !");
 #endif
 
+/* Do we have gravity accuracy checks ? */
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  if (myrank == 0)
+    message(
+        "WARNING: Checking 1/%d of all gpart for gravity accuracy. Code will "
+        "be slower !",
+        SWIFT_GRAVITY_FORCE_CHECKS);
+#endif
+
   /* Do we choke on FP-exceptions ? */
   if (with_fp_exceptions) {
     feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);

src/Makefile.am

@@ -55,7 +55,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     physical_constants.c potential.c hydro_properties.c \
     runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \
     statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
-    part_type.c xmf.c gravity_properties.c
+    part_type.c xmf.c gravity_properties.c gravity.c
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \

src/const.h

@@ -41,7 +41,7 @@
 
 /* Self gravity stuff. */
 #define const_gravity_multipole_order 1
-#define const_gravity_a_smooth 1.25f
+#define const_gravity_a_smooth FLT_MAX
 #define const_gravity_r_cut 4.5f
 #define const_gravity_eta 0.025f
 

src/engine.c

@@ -55,6 +55,7 @@
 #include "cycle.h"
 #include "debug.h"
 #include "error.h"
+#include "gravity.h"
 #include "hydro.h"
 #include "minmax.h"
 #include "parallel_io.h"
@@ -3066,11 +3067,21 @@ void engine_step(struct engine *e) {
   /* 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
+
   /* Start all the tasks. */
   TIMER_TIC;
   engine_launch(e, e->nr_threads);
   TIMER_TOC(timer_runners);
 
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  /* Check the accuracy of the gravity calculation */
+  gravity_exact_force_check(e->s, e, 1e-1);
+#endif
+
 /* Collect the values of rebuild from all nodes. */
 #ifdef WITH_MPI
   int buff = 0;

src/gravity.h

@@ -22,9 +22,20 @@
 /* Config parameters. */
 #include "../config.h"
 
+/* Local headers. */
+#include "const.h"
+#include "engine.h"
+#include "inline.h"
+#include "part.h"
+#include "space.h"
+
 /* So far only one model here */
 /* Straight-forward import */
 #include "./gravity/Default/gravity.h"
 #include "./gravity/Default/gravity_iact.h"
 
+void gravity_exact_force_compute(struct space *s, const struct engine *e);
+void gravity_exact_force_check(struct space *s, const struct engine *e,
+                               float rel_tol);
+
 #endif

src/gravity/Default/gravity_part.h

@@ -63,6 +63,13 @@ struct gpart {
 
 #endif
 
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+
+  /* Brute-force particle acceleration. */
+  float a_grav_exact[3];
+
+#endif
+
 } SWIFT_STRUCT_ALIGN;
 
 #endif /* SWIFT_DEFAULT_GRAVITY_PART_H */