diff --git a/src/debug.c b/src/debug.c
index 24b7938412d05c42b6169939e3f77fe4c755e9af..3732ee5e769277deb393926ea2dc6f04fba93782 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -49,7 +49,7 @@
 #include "./hydro/Default/hydro_debug.h"
 #elif defined(GIZMO_SPH)
 #include "./hydro/Gizmo/hydro_debug.h"
-#elif defined(SHADOWSWIFT)
+#elif defined(SHADOWFAX_SPH)
 #include "./hydro/Shadowswift/hydro_debug.h"
 #else
 #error "Invalid choice of SPH variant"
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 5a37b4144cdffff948d10c717ee6e2fb8b7fda48..25d8b30f8b643c7357286f5df7ea1414648160b7 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -42,6 +42,23 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
   return CFL_condition * R / fabsf(p->timestepvars.vmax);
 }
 
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * We use this to store the physical time step, since it is used for the flux
+ * exchange during the force loop.
+ *
+ * We also set the active flag of the particle to inactive. It will be set to
+ * active in hydro_init_part, which is called the next time the particle becomes
+ * active.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+    struct part* p, float dt) {}
+
 /**
  * @brief Initialises the particles for the first time
  *
@@ -99,10 +116,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
  * This method also initializes the gradient variables (if gradients are used).
  *
  * @param p The particle to act upon.
- * @param The current physical time.
  */
 __attribute__((always_inline)) INLINE static void hydro_end_density(
-    struct part* restrict p, float time) {
+    struct part* restrict p) {
 
   float volume;
   float m, momentum[3], energy;
@@ -164,11 +180,10 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
  * @param timeBase Conversion factor between integer time and physical time.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
-    struct part* restrict p, struct xpart* restrict xp, int ti_current,
-    double timeBase) {
+    struct part* restrict p, struct xpart* restrict xp) {
 
   /* Set the physical time step */
-  p->force.dt = (p->ti_end - p->ti_begin) * timeBase;
+  //  p->force.dt = (p->ti_end - p->ti_begin) * timeBase;
 
   /* Initialize time step criterion variables */
   p->timestepvars.vmax = 0.0f;
@@ -213,6 +228,16 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   p->force.h_dt = 0.0f;
 }
 
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param p The particle.
+ * @param xp The extended data of this particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
+    struct part* restrict p, const struct xpart* restrict xp) {}
+
 /**
  * @brief Converts the hydrodynamic variables from the initial condition file to
  * conserved variables that can be used during the integration
@@ -228,9 +253,10 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
  * velocity of the particle.
  *
  * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
  */
 __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
-    struct part* p) {
+    struct part* p, struct xpart* xp) {
 
   float volume;
   float m;
@@ -258,13 +284,9 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
  * @param p Particle to act upon.
  * @param xp The extended particle data to act upon.
  * @param dt The drift time-step.
- * @param t0 Integer start time of the drift interval.
- * @param t1 Integer end time of the drift interval.
- * @param timeBase Conversion factor between integer and physical time.
  */
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(
-    struct part* p, struct xpart* xp, float dt, int t0, int t1,
-    double timeBase) {}
+    struct part* p, struct xpart* xp, float dt) {}
 
 /**
  * @brief Set the particle acceleration after the flux loop
@@ -361,20 +383,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
  * @param p Particle to act upon.
  * @param xp Extended particle data to act upon.
  * @param dt Physical time step.
- * @param half_dt Half the physical time step.
  */
 __attribute__((always_inline)) INLINE static void hydro_kick_extra(
-    struct part* p, struct xpart* xp, float dt, float half_dt) {}
+    struct part* p, struct xpart* xp, float dt) {}
 
 /**
  * @brief Returns the internal energy of a particle
  *
  * @param p The particle of interest.
- * @param dt Time since the last kick.
  * @return Internal energy of the particle.
  */
 __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
-    const struct part* restrict p, float dt) {
+    const struct part* restrict p) {
 
   return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
 }
@@ -383,11 +403,10 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
  * @brief Returns the entropy of a particle
  *
  * @param p The particle of interest.
- * @param dt Time since the last kick.
  * @return Entropy of the particle.
  */
 __attribute__((always_inline)) INLINE static float hydro_get_entropy(
-    const struct part* restrict p, float dt) {
+    const struct part* restrict p) {
 
   return p->primitives.P / pow_gamma(p->primitives.rho);
 }
@@ -396,11 +415,10 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
  * @brief Returns the sound speed of a particle
  *
  * @param p The particle of interest.
- * @param dt Time since the last kick.
  * @param Sound speed of the particle.
  */
 __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
-    const struct part* restrict p, float dt) {
+    const struct part* restrict p) {
 
   return sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
 }
@@ -409,11 +427,10 @@ __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
  * @brief Returns the pressure of a particle
  *
  * @param p The particle of interest
- * @param dt Time since the last kick
  * @param Pressure of the particle.
  */
 __attribute__((always_inline)) INLINE static float hydro_get_pressure(
-    const struct part* restrict p, float dt) {
+    const struct part* restrict p) {
 
   return p->primitives.P;
 }
diff --git a/src/hydro/Shadowswift/hydro_io.h b/src/hydro/Shadowswift/hydro_io.h
index 3c4037536c3c8dc1012edbf17b3e73998e30b223..c7b1b3fe74ecc8899a7675700ab661a366786968 100644
--- a/src/hydro/Shadowswift/hydro_io.h
+++ b/src/hydro/Shadowswift/hydro_io.h
@@ -143,18 +143,18 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
  */
 void writeSPHflavour(hid_t h_grpsph) {
   /* Gradient information */
-  writeAttribute_s(h_grpsph, "Gradient reconstruction model",
-                   HYDRO_GRADIENT_IMPLEMENTATION);
+  io_write_attribute_s(h_grpsph, "Gradient reconstruction model",
+                       HYDRO_GRADIENT_IMPLEMENTATION);
 
   /* Slope limiter information */
-  writeAttribute_s(h_grpsph, "Cell wide slope limiter model",
-                   HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION);
-  writeAttribute_s(h_grpsph, "Piecewise slope limiter model",
-                   HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION);
+  io_write_attribute_s(h_grpsph, "Cell wide slope limiter model",
+                       HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION);
+  io_write_attribute_s(h_grpsph, "Piecewise slope limiter model",
+                       HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION);
 
   /* Riemann solver information */
-  writeAttribute_s(h_grpsph, "Riemann solver type",
-                   RIEMANN_SOLVER_IMPLEMENTATION);
+  io_write_attribute_s(h_grpsph, "Riemann solver type",
+                       RIEMANN_SOLVER_IMPLEMENTATION);
 }
 
 /**
diff --git a/src/hydro_io.h b/src/hydro_io.h
index b2832bc72e0577fcb0152be5bd91da23bc583512..202d724f821b570b427210bb48b7070563513458 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -32,7 +32,7 @@
 #include "./hydro/Default/hydro_io.h"
 #elif defined(GIZMO_SPH)
 #include "./hydro/Gizmo/hydro_io.h"
-#elif defined(SHADOWSWIFT)
+#elif defined(SHADOWFAX_SPH)
 #include "./hydro/Shadowswift/hydro_io.h"
 #else
 #error "Invalid choice of SPH variant"
diff --git a/tests/benchmarkInteractions.c b/tests/benchmarkInteractions.c
index e3f558f88dffbab252bf7c06f9e943ff568b6fff..0f5b3d2eb294c13e3035885b511c702e6f0cd540 100644
--- a/tests/benchmarkInteractions.c
+++ b/tests/benchmarkInteractions.c
@@ -92,7 +92,7 @@ struct part *make_particles(size_t count, double *offset, double spacing,
   p->h = h;
   p->id = ++(*partId);
 
-#if !defined(GIZMO_SPH)
+#if !defined(GIZMO_SPH) && !defined(SHADOWFAX_SPH)
   p->mass = 1.0f;
 #endif
 
@@ -113,7 +113,7 @@ struct part *make_particles(size_t count, double *offset, double spacing,
 
     p->h = h;
     p->id = ++(*partId);
-#if !defined(GIZMO_SPH)
+#if !defined(GIZMO_SPH) && !defined(SHADOWFAX_SPH)
     p->mass = 1.0f;
 #endif
   }
@@ -125,7 +125,7 @@ struct part *make_particles(size_t count, double *offset, double spacing,
  */
 void prepare_force(struct part *parts, size_t count) {
 
-#if !defined(GIZMO_SPH)
+#if !defined(GIZMO_SPH) && !defined(SHADOWFAX_SPH)
   struct part *p;
   for (size_t i = 0; i < count; ++i) {
     p = &parts[i];
@@ -152,18 +152,18 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
           "%13e %13e %13e\n",
           p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2],
           p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
-#if defined(GIZMO_SPH)
+#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
           0., 0.,
 #else
           p->rho, p->density.rho_dh,
 #endif
           p->density.wcount, p->density.wcount_dh, p->force.h_dt,
-#if defined(GIZMO_SPH)
+#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
           0.,
 #else
           p->force.v_sig,
 #endif
-#if defined(MINIMAL_SPH) || defined(GIZMO_SPH)
+#if defined(MINIMAL_SPH) || defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
           0., 0., 0., 0.
 #else
           p->density.div_v, p->density.rot_v[0], p->density.rot_v[1],
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 44869291d0a69bb22f1dcd50e0fc747368f6843e..aa2b577d606e97b0c2466f420870fa2855fd722e 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -101,7 +101,7 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
   part->u = pressure / (hydro_gamma_minus_one * density);
 #elif defined(MINIMAL_SPH)
   part->u = pressure / (hydro_gamma_minus_one * density);
-#elif defined(GIZMO_SPH) || defined(SHADOWSWIFT)
+#elif defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
   part->primitives.P = pressure;
 #else
   error("Need to define pressure here !");
@@ -198,7 +198,7 @@ void reset_particles(struct cell *c, enum velocity_field vel,
     set_velocity(p, vel, size);
     set_energy_state(p, press, size, density);
 
-#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
+#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
     float volume = p->conserved.mass / density;
 #if defined(GIZMO_SPH)
     p->geometry.volume = volume;
@@ -265,7 +265,7 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
         part->x[2] = offset[2] + size * (z + 0.5) / (float)n;
         part->h = size * h / (float)n;
 
-#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
+#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
         part->conserved.mass = density * volume / count;
 #else
         part->mass = density * volume / count;
@@ -368,7 +368,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
             main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
             main_cell->parts[pid].v[2], main_cell->parts[pid].h,
             hydro_get_density(&main_cell->parts[pid]),
-#if defined(MINIMAL_SPH) || defined(SHADOWSWIFT)
+#if defined(MINIMAL_SPH) || defined(SHADOWFAX_SPH)
             0.f,
 #else
             main_cell->parts[pid].density.div_v,
@@ -422,10 +422,6 @@ void runner_doself1_density(struct runner *r, struct cell *ci);
 void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
 void runner_doself2_force(struct runner *r, struct cell *ci);
 
-#if defined(SHADOWSWIFT) && defined(HYDRO_DIMENSION_3D)
-VORONOI3D_DECLARE_GLOBAL_VARIABLES()
-#endif
-
 /* And go... */
 int main(int argc, char *argv[]) {
 
@@ -446,12 +442,6 @@ int main(int argc, char *argv[]) {
   /* Get some randomness going */
   srand(0);
 
-#if defined(SHADOWSWIFT) && defined(HYDRO_DIMENSION_3D)
-  float box_anchor[3] = {-2.0f, -2.0f, -2.0f};
-  float box_side[3] = {8.0f, 8.0f, 8.0f};
-  voronoi_set_box(box_anchor, box_side);
-#endif
-
   char c;
   while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:p:")) != -1) {
     switch (c) {
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 96a1353be020923282ff93a782e0f52bc873298b..bece9b5e7ac5af547c8f14392ab544e371666869 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -216,7 +216,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
             main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
             main_cell->parts[pid].v[2],
             hydro_get_density(&main_cell->parts[pid]),
-#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
+#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
             0.f,
 #else
             main_cell->parts[pid].density.rho_dh,
@@ -253,7 +253,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
               cj->parts[pjd].id, cj->parts[pjd].x[0], cj->parts[pjd].x[1],
               cj->parts[pjd].x[2], cj->parts[pjd].v[0], cj->parts[pjd].v[1],
               cj->parts[pjd].v[2], hydro_get_density(&cj->parts[pjd]),
-#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
+#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
               0.f,
 #else
               main_cell->parts[pjd].density.rho_dh,
@@ -299,10 +299,6 @@ void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
 void runner_doself1_density(struct runner *r, struct cell *ci);
 void runner_doself1_density_vec(struct runner *r, struct cell *ci);
 
-#if defined(SHADOWSWIFT) && defined(HYDRO_DIMENSION_3D)
-VORONOI3D_DECLARE_GLOBAL_VARIABLES()
-#endif
-
 /* And go... */
 int main(int argc, char *argv[]) {
 
@@ -325,12 +321,6 @@ int main(int argc, char *argv[]) {
   /* Get some randomness going */
   srand(0);
 
-#if defined(SHADOWSWIFT) && defined(HYDRO_DIMENSION_3D)
-  float box_anchor[3] = {-2.0f, -2.0f, -2.0f};
-  float box_side[3] = {6.0f, 6.0f, 6.0f};
-  voronoi_set_box(box_anchor, box_side);
-#endif
-
   char c;
   while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:a:")) != -1) {
     switch (c) {
diff --git a/tests/testPair.c b/tests/testPair.c
index 57916023209a3400bf70831a6ef9e8eaa0a6cca3..6a0a53f7d4a34f4b055080f9c26a394a199b89b4 100644
--- a/tests/testPair.c
+++ b/tests/testPair.c
@@ -142,7 +142,7 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
             ci->parts[pid].id, ci->parts[pid].x[0], ci->parts[pid].x[1],
             ci->parts[pid].x[2], ci->parts[pid].v[0], ci->parts[pid].v[1],
             ci->parts[pid].v[2], hydro_get_density(&ci->parts[pid]),
-#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
+#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
             0.f,
 #else
             ci->parts[pid].density.rho_dh,
@@ -166,7 +166,7 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
             cj->parts[pjd].id, cj->parts[pjd].x[0], cj->parts[pjd].x[1],
             cj->parts[pjd].x[2], cj->parts[pjd].v[0], cj->parts[pjd].v[1],
             cj->parts[pjd].v[2], hydro_get_density(&cj->parts[pjd]),
-#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
+#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
             0.f,
 #else
             cj->parts[pjd].density.rho_dh,
@@ -187,10 +187,6 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
 /* Just a forward declaration... */
 void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
 
-#if defined(SHADOWSWIFT) && defined(HYDRO_DIMENSION_3D)
-VORONOI3D_DECLARE_GLOBAL_VARIABLES()
-#endif
-
 int main(int argc, char *argv[]) {
   size_t particles = 0, runs = 0, volume, type = 0;
   double offset[3] = {0, 0, 0}, h = 1.1255, size = 1., rho = 1.;
@@ -211,12 +207,6 @@ int main(int argc, char *argv[]) {
 
   srand(0);
 
-#if defined(SHADOWSWIFT) && defined(HYDRO_DIMENSION_3D)
-  float box_anchor[3] = {-2.0f, -2.0f, -2.0f};
-  float box_side[3] = {6.0f, 6.0f, 6.0f};
-  voronoi_set_box(box_anchor, box_side);
-#endif
-
   while ((c = getopt(argc, argv, "h:p:r:t:d:f:")) != -1) {
     switch (c) {
       case 'h':