Fixed bug in 2D Voronoi algorithm. Sod shock now works (with fixed cells)....

Fixed bug in 2D Voronoi algorithm. Sod shock now works (with fixed cells). Moved gradient and slope limiter defines for SHADOWFAX_SPH to const.h.

src/const.h

@@ -57,9 +57,17 @@
 //#define GIZMO_FIX_PARTICLES
 //#define GIZMO_TOTAL_ENERGY
 
+/* Types of gradients to use for SHADOWFAX_SPH */
+/* If no option is chosen, no gradients are used (first order scheme) */
+#define SHADOWFAX_GRADIENTS
+
+/* SHADOWFAX_SPH slope limiters */
+#define SHADOWFAX_SLOPE_LIMITER_PER_FACE
+#define SHADOWFAX_SLOPE_LIMITER_CELL_WIDE
+
 /* Options to control SHADOWFAX_SPH */
 /* This option disables cell movement */
-//#define SHADOWFAX_FIX_CELLS
+#define SHADOWFAX_FIX_CELLS
 /* This option evolves the total energy instead of the thermal energy */
 //#define SHADOWFAX_TOTAL_ENERGY
 

src/hydro/Shadowswift/hydro_gradients.h

@@ -20,8 +20,6 @@
 #ifndef SWIFT_HYDRO_GRADIENTS_H
 #define SWIFT_HYDRO_GRADIENTS_H
 
-#define SHADOWFAX_GRADIENTS
-
 #include "hydro_slope_limiters.h"
 
 #if defined(SHADOWFAX_GRADIENTS)

src/hydro/Shadowswift/hydro_slope_limiters.h

@@ -20,13 +20,10 @@
 #ifndef SWIFT_HYDRO_SLOPE_LIMITERS_H
 #define SWIFT_HYDRO_SLOPE_LIMITERS_H
 
-#define PER_FACE_LIMITER
-#define CELL_WIDE_LIMITER
-
 #include "dimension.h"
 #include "kernel_hydro.h"
 
-#ifdef PER_FACE_LIMITER
+#ifdef SHADOWFAX_SLOPE_LIMITER_PER_FACE
 
 #define HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION \
   "GIZMO piecewise slope limiter (Hopkins 2015)"
@@ -55,7 +52,7 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
 
 #endif
 
-#ifdef CELL_WIDE_LIMITER
+#ifdef SHADOWFAX_SLOPE_LIMITER_CELL_WIDE
 
 #define HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION \
   "Cell wide slope limiter (Springel 2010)"

src/hydro/Shadowswift/voronoi2d_algorithm.h

@@ -479,8 +479,9 @@ __attribute__((always_inline)) INLINE float voronoi_cell_finalize(
     cell->centroid[0] += (x[0] + x[1]) * A;
     cell->centroid[1] += (y[0] + y[1]) * A;
 
-    cell->face_midpoints[i][0] = 0.5f * (x[0] + x[1]) + cell->x[0];
-    cell->face_midpoints[i][1] = 0.5f * (y[0] + y[1]) + cell->x[1];
+    /* Note that we only need the RELATIVE positions of the midpoints */
+    cell->face_midpoints[i][0] = 0.5f * (x[0] + x[1]);
+    cell->face_midpoints[i][1] = 0.5f * (y[0] + y[1]);
 
     r2 = x[0] * x[0] + y[0] * y[0];
     r2max = max(r2max, r2);