Make the code crash if the wrong softening kernel is used in the high-order derivatives

src/gravity_properties.c

@@ -114,6 +114,13 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
   p->use_tree_below_softening =
       parser_get_opt_param_int(params, "Gravity:use_tree_below_softening", 1);
 
+#ifdef GADGET2_SOFTENING_CORRECTION
+  if (p->use_tree_below_softening)
+    error(
+        "Cannot solve gravity via the tree below softening with the "
+        "Gadget2-type softening kernel");
+#endif
+
   /* Mesh dithering */
   if (periodic && !with_external_potential) {
     p->with_dithering =

src/kernel_gravity.h

@@ -162,8 +162,6 @@ __attribute__((always_inline)) INLINE static void kernel_grav_eval_force_double(
 }
 #endif /* SWIFT_GRAVITY_FORCE_CHECKS */
 
-#undef GADGET2_SOFTENING_CORRECTION
-
 /************************************************/
 /* Derivatives of softening kernel used for FMM */
 /************************************************/
@@ -171,6 +169,10 @@ __attribute__((always_inline)) INLINE static void kernel_grav_eval_force_double(
 __attribute__((const)) INLINE static float D_soft_1(const float u,
                                                     const float u_inv) {
 
+#ifdef GADGET2_SOFTENING_CORRECTION
+  error("Invalid choice of softening kernel shape");
+#endif
+
   /* phi(u) = 3u^7 - 15u^6 + 28u^5 - 21u^4 + 7u^2 - 3 */
   float phi = 3.f * u - 15.f;
   phi = phi * u + 28.f;
@@ -186,6 +188,10 @@ __attribute__((const)) INLINE static float D_soft_1(const float u,
 __attribute__((const)) INLINE static float D_soft_3(const float u,
                                                     const float u_inv) {
 
+#ifdef GADGET2_SOFTENING_CORRECTION
+  error("Invalid choice of softening kernel shape");
+#endif
+
   /* phi'(u)/u = 21u^5 - 90u^4 + 140u^3 - 84u^2 + 14 */
   float phi = 21.f * u - 90.f;
   phi = phi * u + 140.f;
@@ -199,6 +205,10 @@ __attribute__((const)) INLINE static float D_soft_3(const float u,
 __attribute__((const)) INLINE static float D_soft_5(const float u,
                                                     const float u_inv) {
 
+#ifdef GADGET2_SOFTENING_CORRECTION
+  error("Invalid choice of softening kernel shape");
+#endif
+
   /* (phi'(u)/u)'/u = 105u^3 - 360u^2 + 420u - 168 */
   float phi = 105.f * u - 360.f;
   phi = phi * u + 420.f;
@@ -209,11 +219,20 @@ __attribute__((const)) INLINE static float D_soft_5(const float u,
 
 __attribute__((const)) INLINE static float D_soft_7(const float u,
                                                     const float u_inv) {
+#ifdef GADGET2_SOFTENING_CORRECTION
+  error("Invalid choice of softening kernel shape");
+#endif
+
   return 315.f * u - 720.f + 420.f * u_inv;
 }
 
 __attribute__((const)) INLINE static float D_soft_9(const float u,
                                                     const float u_inv) {
+
+#ifdef GADGET2_SOFTENING_CORRECTION
+  error("Invalid choice of softening kernel shape");
+#endif
+
   /* -315 u^-1 + 420 u^-3 */
   /* -315 u_inv + 420 u_inv^3 */
   float phi = 420.f * u_inv;
@@ -224,6 +243,11 @@ __attribute__((const)) INLINE static float D_soft_9(const float u,
 
 __attribute__((const)) INLINE static float D_soft_11(const float u,
                                                      const float u_inv) {
+
+#ifdef GADGET2_SOFTENING_CORRECTION
+  error("Invalid choice of softening kernel shape");
+#endif
+
   /* 315 u^-3 - 1260 u^-5 */
   /* 315 u_inv^3 - 1260 u_inv^5 */
   float phi = -1260.f * u_inv;