diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index e4ea0d971cd89862d3ea6b99f0b04930666a154a..36078798cdd1ac68a456134fd7887408752f18c9 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -20,6 +20,7 @@
 #include <float.h>
 #include "adiabatic_index.h"
 #include "approx_math.h"
+#include "cosmology.h"
 #include "equation_of_state.h"
 #include "hydro_gradients.h"
 #include "hydro_properties.h"
@@ -35,7 +36,8 @@
  */
 __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
     const struct part* restrict p, const struct xpart* restrict xp,
-    const struct hydro_props* restrict hydro_properties) {
+    const struct hydro_props* restrict hydro_properties,
+    const struct cosmology* restrict cosmo) {
 
   const float CFL_condition = hydro_properties->CFL_condition;
 
@@ -158,7 +160,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
  * @param p The particle to act upon.
  */
 __attribute__((always_inline)) INLINE static void hydro_end_density(
-    struct part* restrict p) {
+    struct part* restrict p, const struct cosmology* cosmo) {
 
   float volume;
   float m, momentum[3], energy;
@@ -246,7 +248,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
  * @param xp The extended particle data to act upon
  */
 __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
-    struct part* restrict p, struct xpart* restrict xp) {
+    struct part* restrict p, struct xpart* restrict xp,
+    const struct cosmology* cosmo) {
 
   /* Some smoothing length multiples. */
   const float h = p->h;
@@ -273,7 +276,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @param xp The extended particle data to act upon.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
-    struct part* restrict p, struct xpart* restrict xp) {
+    struct part* restrict p, struct xpart* restrict xp,
+    const struct cosmology* cosmo) {
 
   /* Initialize time step criterion variables */
   p->timestepvars.vmax = 0.0f;
@@ -346,7 +350,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
  * @param xp The extended particle data to act upon.
  */
 __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
-    struct part* p, struct xpart* xp) {}
+    struct part* p, struct xpart* xp, const struct cosmology* cosmo) {}
 
 /**
  * @brief Extra operations to be done during the drift
@@ -358,7 +362,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
  * @param dt The drift time-step.
  */
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(
-    struct part* p, struct xpart* xp, float dt) {}
+    struct part* p, struct xpart* xp, float dt_drift, float dt_therm) {}
 
 /**
  * @brief Set the particle acceleration after the flux loop.
@@ -366,7 +370,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
  * @param p Particle to act upon.
  */
 __attribute__((always_inline)) INLINE static void hydro_end_force(
-    struct part* p) {}
+    struct part* p, const struct cosmology* cosmo) {}
 
 /**
  * @brief Extra operations done during the kick
@@ -378,7 +382,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
  * @param dt Physical time step.
  */
 __attribute__((always_inline)) INLINE static void hydro_kick_extra(
-    struct part* p, struct xpart* xp, float dt) {
+    struct part* p, struct xpart* xp, float dt, const struct cosmology* cosmo,
+    const struct hydro_props* hydro_props) {
 
   float vcell[3];
 
@@ -553,8 +558,13 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
  * @param v (return) The velocities at the current time.
  */
 __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
-    const struct part* restrict p, const struct xpart* xp, float dt,
-    float v[3]) {}
+    const struct part* restrict p, const struct xpart* xp, float dt_kick_hydro,
+    float dt_kick_grav, float v[3]) {
+
+  v[0] = p->v[0];
+  v[1] = p->v[1];
+  v[2] = p->v[2];
+}
 
 /**
  * @brief Returns the density of a particle
@@ -621,3 +631,171 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
     p->primitives.P = gas_pressure_from_entropy(p->primitives.rho, S);
   }
 }
+
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+    struct part* restrict p, float m) {
+
+  p->conserved.mass = m;
+}
+
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part* p, float u_init) {
+
+  p->conserved.energy = u_init * p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+  /* add the kinetic energy */
+  p->conserved.energy += 0.5f * p->conserved.mass *
+                         (p->conserved.momentum[0] * p->primitives.v[0] +
+                          p->conserved.momentum[1] * p->primitives.v[1] +
+                          p->conserved.momentum[2] * p->primitives.v[2]);
+#endif
+  p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u_init;
+}
+
+/**
+ * @brief Returns the comoving internal energy of a particle
+ *
+ * @param p The particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy(const struct part* restrict p) {
+
+  if (p->primitives.rho > 0.)
+    return gas_internal_energy_from_pressure(p->primitives.rho,
+                                             p->primitives.P);
+  else
+    return 0.;
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle
+ *
+ * @param p The particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
+    const struct part* restrict p) {
+
+  if (p->primitives.rho > 0.) {
+    return gas_entropy_from_pressure(p->primitives.rho, p->primitives.P);
+  } else {
+    return 0.;
+  }
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_soundspeed(const struct part* restrict p) {
+
+  if (p->primitives.rho > 0.)
+    return gas_soundspeed_from_pressure(p->primitives.rho, p->primitives.P);
+  else
+    return 0.;
+}
+
+/**
+ * @brief Returns the comoving pressure of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
+    const struct part* restrict p) {
+
+  return p->primitives.P;
+}
+
+/**
+ * @brief Returns the comoving density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
+    const struct part* restrict p) {
+
+  return p->primitives.rho;
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy(const struct part* restrict p,
+                                   const struct cosmology* cosmo) {
+
+  return cosmo->a_factor_internal_energy *
+         hydro_get_comoving_internal_energy(p);
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
+    const struct part* restrict p, const struct cosmology* cosmo) {
+
+  /* Note: no cosmological conversion required here with our choice of
+   * coordinates. */
+  return hydro_get_comoving_entropy(p);
+}
+
+/**
+ * @brief Returns the physical sound speed of a particle
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_soundspeed(const struct part* restrict p,
+                              const struct cosmology* cosmo) {
+
+  return cosmo->a_factor_sound_speed * hydro_get_comoving_soundspeed(p);
+}
+
+/**
+ * @brief Returns the comoving pressure of a particle
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
+    const struct part* restrict p, const struct cosmology* cosmo) {
+
+  return cosmo->a_factor_pressure * p->primitives.P;
+}
+
+/**
+ * @brief Returns the physical density of a particle
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_density(
+    const struct part* restrict p, const struct cosmology* cosmo) {
+
+  return cosmo->a3_inv * p->primitives.rho;
+}
diff --git a/src/hydro/Shadowswift/hydro_gradients.h b/src/hydro/Shadowswift/hydro_gradients.h
index 1aea49790d998d3912a80fa1376cbd1e183f26f7..285d889a1a6e10662a06979f69290aabd4206059 100644
--- a/src/hydro/Shadowswift/hydro_gradients.h
+++ b/src/hydro/Shadowswift/hydro_gradients.h
@@ -51,8 +51,8 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_init(
  * @param pj Particle j.
  */
 __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
-    float r2, float* dx, float hi, float hj, struct part* pi, struct part* pj) {
-}
+    float r2, const float* dx, float hi, float hj, struct part* restrict pi,
+    struct part* restrict pj) {}
 
 /**
  * @brief Gradient calculations done during the neighbour loop: non-symmetric
@@ -66,8 +66,9 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
  * @param pj Particle j.
  */
 __attribute__((always_inline)) INLINE static void
-hydro_gradients_nonsym_collect(float r2, float* dx, float hi, float hj,
-                               struct part* pi, struct part* pj) {}
+hydro_gradients_nonsym_collect(float r2, const float* dx, float hi, float hj,
+                               struct part* restrict pi,
+                               const struct part* restrict pj) {}
 
 /**
  * @brief Finalize the gradient variables after all data have been collected
@@ -84,8 +85,8 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
  * gradients_none does nothing, since all gradients are zero -- are they?).
  */
 __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
-    struct part* pi, struct part* pj, float hi, float hj, float* dx, float r,
-    float* xij_i, float* Wi, float* Wj, float mindt) {
+    struct part* pi, struct part* pj, float hi, float hj, const float* dx,
+    float r, float* xij_i, float* Wi, float* Wj, float mindt) {
 
   float dWi[5], dWj[5];
   float xij_j[3];
diff --git a/src/hydro/Shadowswift/hydro_gradients_shadowfax.h b/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
index 9ca40a604da3dc12bbb48ac033cd078f0561d8ab..d131731907806536e86a03921b1c701f287077f1 100644
--- a/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
+++ b/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
@@ -65,7 +65,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_init(
  * @param grad Current value of the gradient for the quantity (is updated).
  */
 __attribute__((always_inline)) INLINE void hydro_gradients_single_quantity(
-    float qL, float qR, float *cLR, float *xLR, float rLR, float A,
+    float qL, float qR, float *cLR, const float *xLR, float rLR, float A,
     float *grad) {
 
   grad[0] += A * ((qR - qL) * cLR[0] / rLR - 0.5f * (qL + qR) * xLR[0] / rLR);
@@ -84,7 +84,8 @@ __attribute__((always_inline)) INLINE void hydro_gradients_single_quantity(
  * @param pj Particle j.
  */
 __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+    float r2, const float *dx, float hi, float hj, struct part *pi,
+    struct part *pj) {
 
   float A, midpoint[3];
 
@@ -148,8 +149,8 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
  * @param pj Particle j.
  */
 __attribute__((always_inline)) INLINE static void
-hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
-                               struct part *pi, struct part *pj) {
+hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
+                               struct part *pi, const struct part *pj) {
 
   float A, midpoint[3];
 
diff --git a/src/hydro/Shadowswift/hydro_iact.h b/src/hydro/Shadowswift/hydro_iact.h
index 15219bf879ea7d34a8a45be217ef81107d3392e6..9ac1debf3184c25603412867c41c62a1131345f3 100644
--- a/src/hydro/Shadowswift/hydro_iact.h
+++ b/src/hydro/Shadowswift/hydro_iact.h
@@ -35,7 +35,8 @@
  * @param pj Particle j.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_density(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {
 
   float mindx[3];
 
@@ -59,7 +60,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
  * @param pj Particle j.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    const struct part *restrict pj, float a, float H) {
 
   voronoi_cell_interact(&pi->cell, dx, pj->id);
 }
@@ -78,7 +80,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
  * @param pj Particle j.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_gradient(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {
 
   hydro_gradients_collect(r2, dx, hi, hj, pi, pj);
 }
@@ -98,7 +101,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
  * @param pj Particle j.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {
 
   hydro_gradients_nonsym_collect(r2, dx, hi, hj, pi, pj);
 }
@@ -129,8 +133,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
  * @param pj Particle j.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
-    int mode) {
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, int mode, float a, float H) {
 
   float r = sqrtf(r2);
   int k;
@@ -322,9 +326,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
  * @param pj Particle j.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_force(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {
 
-  runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1);
+  runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1, a, H);
 }
 
 /**
@@ -341,7 +346,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
  * @param pj Particle j.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {
 
-  runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0);
+  runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0, a, H);
 }
diff --git a/src/hydro/Shadowswift/hydro_io.h b/src/hydro/Shadowswift/hydro_io.h
index 65cb4ee6b8b1ecb174212bc8867cd7657a213414..8525d22025c1943529ddcd86cf3a42ba0ae4f5d4 100644
--- a/src/hydro/Shadowswift/hydro_io.h
+++ b/src/hydro/Shadowswift/hydro_io.h
@@ -64,7 +64,8 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
  * @param p Particle.
  * @return Internal energy of the particle
  */
-void convert_u(const struct engine* e, const struct part* p, float* ret) {
+void convert_u(const struct engine* e, const struct part* p,
+               const struct xpart* xp, float* ret) {
   ret[0] = hydro_get_internal_energy(p);
 }
 
@@ -75,7 +76,8 @@ void convert_u(const struct engine* e, const struct part* p, float* ret) {
  * @param p Particle.
  * @return Entropic function of the particle
  */
-void convert_A(const struct engine* e, const struct part* p, float* ret) {
+void convert_A(const struct engine* e, const struct part* p,
+               const struct xpart* xp, float* ret) {
   ret[0] = hydro_get_entropy(p);
 }
 
@@ -86,7 +88,8 @@ void convert_A(const struct engine* e, const struct part* p, float* ret) {
  * @param p Particle.
  * @return Total energy of the particle
  */
-void convert_Etot(const struct engine* e, const struct part* p, float* ret) {
+void convert_Etot(const struct engine* e, const struct part* p,
+                  const struct xpart* xp, float* ret) {
 #ifdef SHADOWFAX_TOTAL_ENERGY
   return p->conserved.energy;
 #else
@@ -105,7 +108,7 @@ void convert_Etot(const struct engine* e, const struct part* p, float* ret) {
 }
 
 void convert_part_pos(const struct engine* e, const struct part* p,
-                      double* ret) {
+                      const struct xpart* xp, double* ret) {
 
   if (e->s->periodic) {
     ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
@@ -125,14 +128,15 @@ void convert_part_pos(const struct engine* e, const struct part* p,
  * @param list The list of i/o properties to write.
  * @param num_fields The number of i/o fields to write.
  */
-void hydro_write_particles(struct part* parts, struct io_props* list,
-                           int* num_fields) {
+void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
+                           struct io_props* list, int* num_fields) {
 
   *num_fields = 13;
 
   /* List what we want to write */
-  list[0] = io_make_output_field_convert_part(
-      "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
+  list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+                                              UNIT_CONV_LENGTH, parts, xparts,
+                                              convert_part_pos);
   list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts,
                                  primitives.v);
   list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts,
@@ -141,7 +145,7 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
                                  parts, h);
   list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
                                               UNIT_CONV_ENERGY_PER_UNIT_MASS,
-                                              parts, convert_u);
+                                              parts, xparts, convert_u);
   list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
                                  UNIT_CONV_NO_UNITS, parts, id);
   list[6] = io_make_output_field("Acceleration", FLOAT, 3,
@@ -153,11 +157,11 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
   list[9] = io_make_output_field("GradDensity", FLOAT, 3, UNIT_CONV_DENSITY,
                                  parts, primitives.gradients.rho);
   list[10] = io_make_output_field_convert_part(
-      "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, convert_A);
+      "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, xparts, convert_A);
   list[11] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
                                   parts, primitives.P);
   list[12] = io_make_output_field_convert_part(
-      "TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, convert_Etot);
+      "TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, xparts, convert_Etot);
 }
 
 /**
diff --git a/src/hydro/Shadowswift/hydro_part.h b/src/hydro/Shadowswift/hydro_part.h
index 4f28d7bdc1c4779a9e6f15c3fb789c3f65f32890..a7cc9daf0839216f098ac05c2267adc60ea11fb0 100644
--- a/src/hydro/Shadowswift/hydro_part.h
+++ b/src/hydro/Shadowswift/hydro_part.h
@@ -35,6 +35,9 @@ struct xpart {
   /* Velocity at the last full step. */
   float v_full[3];
 
+  /* Gravitational acceleration at the last full step. */
+  float a_grav[3];
+
   /* Additional data used to record cooling information */
   struct cooling_xpart_data cooling_data;
 
diff --git a/src/hydro/Shadowswift/hydro_slope_limiters_cell.h b/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
index a7b3f7511fa7cd247fd9d2399cd200d8d943630e..d746ffc59538fcc625a2e095245adfd7a946a95e 100644
--- a/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
+++ b/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
@@ -50,7 +50,8 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
  * @param r Distance between particle i and particle j.
  */
 __attribute__((always_inline)) INLINE static void
-hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
+hydro_slope_limit_cell_collect(struct part* pi, const struct part* pj,
+                               float r) {
 
   /* basic slope limiter: collect the maximal and the minimal value for the
    * primitive variables among the ngbs */