Added functions for periodic box wrapping

src/Makefile.am

@@ -63,7 +63,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
 		 kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h runner_doiact_fft.h \
                  runner_doiact_nosort.h units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
-		 dimension.h equation_of_state.h part_type.h \
+		 dimension.h equation_of_state.h part_type.h periodic.h
 		 gravity.h gravity_io.h \
 		 gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
 		 gravity/Default/gravity_debug.h gravity/Default/gravity_part.h  \

src/minmax.h

@@ -43,18 +43,4 @@
     _a > _b ? _a : _b;            \
   })
 
-/**
- * @brief Limits the value of x to be between a and b
- *
- * Only wraps once. If x > 2b, the returned value will be larger than b.
- * Similarly for x < -b.
- */
-#define box_wrap(x, a, b)                               \
-  ({                                                    \
-    const __typeof__(x) _x = (x);                       \
-    const __typeof__(a) _a = (a);                       \
-    const __typeof__(b) _b = (b);                       \
-    _x < _a ? (_x + _b) : ((_x > _b) ? (_x - _b) : _x); \
-  })
-
 #endif /* SWIFT_MINMAX_H */

src/multipole.h

@@ -37,8 +37,8 @@
 #include "gravity_softened_derivatives.h"
 #include "inline.h"
 #include "kernel_gravity.h"
-#include "minmax.h"
 #include "part.h"
+#include "periodic.h"
 #include "vector_power.h"
 
 #define multipole_align 128

src/periodic.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016   Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PERIODIC_H
+#define SWIFT_PERIODIC_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Includes. */
+#include "inline.h"
+
+/**
+ * @brief Limits the value of x to be between a and b
+ *
+ * Only wraps once. If x > 2b, the returned value will be larger than b.
+ * Similarly for x < -b.
+ */
+#define box_wrap(x, a, b)                               \
+  ({                                                    \
+    const __typeof__(x) _x = (x);                       \
+    const __typeof__(a) _a = (a);                       \
+    const __typeof__(b) _b = (b);                       \
+    _x < _a ? (_x + _b) : ((_x > _b) ? (_x - _b) : _x); \
+  })
+
+/**
+ * @brief Find the smallest distance dx along one axis within a box of size
+ * box_size
+ *
+ * This macro evaluates its arguments exactly once.
+ *
+ * Only wraps once. If dx > 2b, the returned value will be larger than b.
+ * Similarly for dx < -b.
+ *
+ */
+__attribute__((always_inline)) INLINE static double nearest(double dx,
+                                                            double box_size) {
+  return dx > 0.5 * box_size ? (dx - box_size)
+                             : ((dx < -0.5 * box_size) ? (dx + box_size) : dx);
+}
+
+/**
+ * @brief Find the smallest distance dx along one axis within a box of size
+ * box_size
+ *
+ * This macro evaluates its arguments exactly once.
+ *
+ * Only wraps once. If dx > 2b, the returned value will be larger than b.
+ * Similarly for dx < -b.
+ *
+ */
+__attribute__((always_inline)) INLINE static float nearestf(float dx,
+                                                            float box_size) {
+  return dx > 0.5f * box_size
+             ? (dx - box_size)
+             : ((dx < -0.5f * box_size) ? (dx + box_size) : dx);
+}
+
+#endif /* SWIFT_PERIODIC_H */

src/swift.h

@@ -45,6 +45,7 @@
 #include "parser.h"
 #include "part.h"
 #include "partition.h"
+#include "periodic.h"
 #include "physical_constants.h"
 #include "potential.h"
 #include "profiler.h"

tests/testFFT.c

@@ -35,8 +35,6 @@ int main() { return 0; }
 /* Includes. */
 #include "swift.h"
 
-#define NEAREST(x) (((x) > 0.5) ? ((x)-1.) : (((x) < -0.5) ? ((x) + 1.) : (x)))
-
 int main() {
 
   /* Initialize CPU frequency, this also starts time. */
@@ -97,17 +95,20 @@ int main() {
   double *pot = malloc(nr_cells * sizeof(double));
   double *pot_exact = malloc(nr_cells * sizeof(double));
 
-  FILE *file = fopen("potential.dat", "w");
+  // FILE *file = fopen("potential.dat", "w");
   for (int i = 0; i < nr_cells; ++i) {
     pot[i] = space.multipoles_top[i].pot.F_000;
-    double dx = NEAREST(space.multipoles_top[i].CoM[0] - gparts[0].x[0]);
-    double dy = NEAREST(space.multipoles_top[i].CoM[1] - gparts[0].x[1]);
-    double dz = NEAREST(space.multipoles_top[i].CoM[2] - gparts[0].x[2]);
+    double dx =
+        nearest(space.multipoles_top[i].CoM[0] - gparts[0].x[0], dim[0]);
+    double dy =
+        nearest(space.multipoles_top[i].CoM[1] - gparts[0].x[1], dim[1]);
+    double dz =
+        nearest(space.multipoles_top[i].CoM[2] - gparts[0].x[2], dim[2]);
     r[i] = sqrt(dx * dx + dy * dy + dz * dz);
     if (r[i] > 0) pot_exact[i] = -1. / r[i];
-    fprintf(file, "%e %e %e\n", r[i], pot[i], pot_exact[i]);
+    // fprintf(file, "%e %e %e\n", r[i], pot[i], pot_exact[i]);
   }
-  fclose(file);
+  // fclose(file);
 
   /* Clean up */
   free(r);