diff --git a/src/Makefile.am b/src/Makefile.am
index 8fca05cd3b8b6a569e0fe75e11a0328b488cbedf..6285d93176c488c4ccf4b948eb886936ffecffc4 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -147,19 +147,19 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  hydro/Gizmo/hydro_unphysical.h \
                  hydro/Gizmo/hydro_gradients_sph.h \
                  hydro/Gizmo/hydro_gradients_gizmo.h \
+                 hydro/Gizmo/hydro_velocities.h \
                  hydro/Gizmo/MFV/hydro_debug.h \
                  hydro/Gizmo/MFV/hydro_part.h \
-                 hydro/Gizmo/MFV/hydro_slope_limiters_cell.h \
                  hydro/Gizmo/MFV/hydro_velocities.h \
                  hydro/Gizmo/MFV/hydro_getters.h \
                  hydro/Gizmo/MFV/hydro_setters.h \
                  hydro/Gizmo/MFV/hydro_flux.h \
                  hydro/Gizmo/MFM/hydro_debug.h \
                  hydro/Gizmo/MFM/hydro_part.h \
-                 hydro/Gizmo/MFM/hydro_slope_limiters_cell.h \
                  hydro/Gizmo/MFM/hydro_getters.h \
                  hydro/Gizmo/MFM/hydro_setters.h \
                  hydro/Gizmo/MFM/hydro_flux.h \
+                 hydro/Gizmo/MFM/hydro_velocities.h \
                  hydro/Shadowswift/hydro_debug.h \
                  hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
                  hydro/Shadowswift/hydro_iact.h \
diff --git a/src/hydro/Gizmo/MFM/hydro_getters.h b/src/hydro/Gizmo/MFM/hydro_getters.h
index 9be5ae04dc92dbcb16c8a4b7e22d364b36cc8cae..e0c348ba9dc5760df674ce859b35c29808f9a25f 100644
--- a/src/hydro/Gizmo/MFM/hydro_getters.h
+++ b/src/hydro/Gizmo/MFM/hydro_getters.h
@@ -19,120 +19,6 @@
 #ifndef SWIFT_GIZMO_MFM_HYDRO_GETTERS_H
 #define SWIFT_GIZMO_MFM_HYDRO_GETTERS_H
 
-#include "cosmology.h"
-#include "equation_of_state.h"
-
-/**
- * @brief Get a 5-element state vector W containing the primitive hydrodynamic
- * variables.
- *
- * @param p Particle.
- * @param W Pointer to the array in which the result needs to be stored (of size
- * 5 or more).
- */
-__attribute__((always_inline)) INLINE static void
-hydro_part_get_primitive_variables(const struct part* restrict p, float* W) {
-
-  W[0] = p->rho;
-  W[1] = p->v[0];
-  W[2] = p->v[1];
-  W[3] = p->v[2];
-  W[4] = p->P;
-}
-
-/**
- * @brief Get the degenerate case correction factor for the given particle.
- *
- * @param p Particle.
- * @return Degenerate case correction factor.
- */
-__attribute__((always_inline)) INLINE static float hydro_part_get_wcorr(
-    const struct part* restrict p) {
-  return p->geometry.wcorr;
-}
-
-/**
- * @brief Get the gradients of the primitive variables for the given particle.
- *
- * @param p Particle.
- * @param drho Density gradient (of size 3 or more).
- * @param dvx x velocity gradient (of size 3 or more).
- * @param dvy y velocity gradient (of size 3 or more).
- * @param dvz z velocity gradient (of size 3 or more).
- * @param dP Pressure gradient (of size 3 or more).
- */
-__attribute__((always_inline)) INLINE static void hydro_part_get_gradients(
-    const struct part* restrict p, float* drho, float* dvx, float* dvy,
-    float* dvz, float* dP) {
-
-  drho[0] = p->gradients.rho[0];
-  drho[1] = p->gradients.rho[1];
-  drho[2] = p->gradients.rho[2];
-
-  dvx[0] = p->gradients.v[0][0];
-  dvx[1] = p->gradients.v[0][1];
-  dvx[2] = p->gradients.v[0][2];
-  dvy[0] = p->gradients.v[1][0];
-  dvy[1] = p->gradients.v[1][1];
-  dvy[2] = p->gradients.v[1][2];
-  dvz[0] = p->gradients.v[2][0];
-  dvz[1] = p->gradients.v[2][1];
-  dvz[2] = p->gradients.v[2][2];
-
-  dP[0] = p->gradients.P[0];
-  dP[1] = p->gradients.P[1];
-  dP[2] = p->gradients.P[2];
-}
-
-/**
- * @brief Get the pressure gradient for the given particle.
- *
- * @param p Particle.
- * @param gradP Pressure gradient (of size 3 or more).
- */
-__attribute__((always_inline)) INLINE static void
-hydro_part_get_pressure_gradient(const struct part* restrict p,
-                                 float gradP[3]) {
-
-  gradP[0] = p->gradients.P[0];
-  gradP[1] = p->gradients.P[1];
-  gradP[2] = p->gradients.P[2];
-}
-
-/**
- * @brief Get the slope limiter variables for the given particle.
- *
- * @param p Particle.
- * @param rholim Minimum and maximum density of neighbours (of size 2 or more).
- * @param vxlim Minimum and maximum x velocity of neighbours (of size 2 or
- * more).
- * @param vylim Minimum and maximum y velocity of neighbours (of size 2 or
- * more).
- * @param vzlim Minimum and maximum z velocity of neighbours (of size 2 or
- * more).
- * @param Plim Minimum and maximum pressure of neighbours (of size 2 or more).
- * @param rmax Maximum distance of any neighbour (of size 1 or more).
- */
-__attribute__((always_inline)) INLINE static void hydro_part_get_slope_limiter(
-    const struct part* restrict p, float* rholim, float* vxlim, float* vylim,
-    float* vzlim, float* Plim, float* rmax) {
-
-  rholim[0] = p->limiter.rho[0];
-  rholim[1] = p->limiter.rho[1];
-
-  vxlim[0] = p->limiter.v[0][0];
-  vxlim[1] = p->limiter.v[0][1];
-  vylim[0] = p->limiter.v[1][0];
-  vylim[1] = p->limiter.v[1][1];
-  vzlim[0] = p->limiter.v[2][0];
-  vzlim[1] = p->limiter.v[2][1];
-
-  Plim[0] = p->limiter.P[0];
-  Plim[1] = p->limiter.P[1];
-
-  rmax[0] = p->limiter.maxr;
-}
-
 /**
  * @brief Get the fluxes for the given particle.
  *
@@ -148,245 +34,4 @@ __attribute__((always_inline)) INLINE static void hydro_part_get_fluxes(
   flux[4] = p->flux.energy;
 }
 
-/**
- * @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->rho > 0.0f) {
-    return gas_internal_energy_from_pressure(p->rho, p->P);
-  } else {
-    return 0.0f;
-  }
-}
-
-/**
- * @brief Returns the physical internal energy of a particle
- *
- * @param p The particle of interest.
- * @param xp The extended data of 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 xpart* restrict xp,
-                                   const struct cosmology* cosmo) {
-
-  return cosmo->a_factor_internal_energy *
-         hydro_get_comoving_internal_energy(p);
-}
-
-/**
- * @brief Returns the comoving internal energy of a particle drifted to the
- * current time.
- *
- * @param p The particle of interest
- */
-__attribute__((always_inline)) INLINE static float
-hydro_get_drifted_comoving_internal_energy(const struct part* restrict p) {
-
-  return hydro_get_comoving_internal_energy(p);
-}
-
-/**
- * @brief Returns the physical internal energy of a particle drifted to the
- * current time.
- *
- * @param p The particle of interest.
- * @param cosmo The cosmological model.
- */
-__attribute__((always_inline)) INLINE static float
-hydro_get_drifted_physical_internal_energy(const struct part* restrict p,
-                                           const struct cosmology* cosmo) {
-
-  return hydro_get_comoving_internal_energy(p) *
-         cosmo->a_factor_internal_energy;
-}
-
-/**
- * @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->rho > 0.0f) {
-    return gas_entropy_from_pressure(p->rho, p->P);
-  } else {
-    return 0.0f;
-  }
-}
-
-/**
- * @brief Returns the physical internal energy of a particle
- *
- * @param p The particle of interest.
- * @param xp The extended data of 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 xpart* restrict xp,
-    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 internal energy of a particle
- *
- * @param p The particle of interest.
- * @param cosmo The cosmological model.
- */
-__attribute__((always_inline)) INLINE static float
-hydro_get_drifted_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 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->rho > 0.0f) {
-    return gas_soundspeed_from_pressure(p->rho, p->P);
-  } else {
-    return 0.0f;
-  }
-}
-
-/**
- * @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
- */
-__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
-    const struct part* restrict p) {
-
-  return p->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->P;
-}
-
-/**
- * @brief Returns the mass of a particle
- *
- * @param p The particle of interest
- */
-__attribute__((always_inline)) INLINE static float hydro_get_mass(
-    const struct part* restrict p) {
-
-  return p->conserved.mass;
-}
-
-/**
- * @brief Returns the velocities drifted to the current time of a particle.
- *
- * @param p The particle of interest
- * @param xp The extended data of the particle.
- * @param dt_kick_hydro The time (for hydro accelerations) since the last kick.
- * @param dt_kick_grav The time (for gravity accelerations) since the last kick.
- * @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_kick_hydro,
-    float dt_kick_grav, float v[3]) {
-
-  if (p->conserved.mass > 0.0f) {
-    const float inverse_mass = 1.0f / p->conserved.mass;
-    v[0] = p->v[0] + p->flux.momentum[0] * dt_kick_hydro * inverse_mass;
-    v[1] = p->v[1] + p->flux.momentum[1] * dt_kick_hydro * inverse_mass;
-    v[2] = p->v[2] + p->flux.momentum[2] * dt_kick_hydro * inverse_mass;
-  } else {
-    v[0] = p->v[0];
-    v[1] = p->v[1];
-    v[2] = p->v[2];
-  }
-
-  v[0] += xp->a_grav[0] * dt_kick_grav;
-  v[1] += xp->a_grav[1] * dt_kick_grav;
-  v[2] += xp->a_grav[2] * dt_kick_grav;
-}
-
-/**
- * @brief Returns the time derivative of co-moving internal energy of a particle
- *
- * We assume a constant density.
- *
- * @param p The particle of interest
- */
-__attribute__((always_inline)) INLINE static float
-hydro_get_comoving_internal_energy_dt(const struct part* restrict p) {
-
-  error("Needs implementing");
-  return 0.f;
-}
-
-/**
- * @brief Returns the time derivative of physical internal energy of a particle
- *
- * We assume a constant density.
- *
- * @param p The particle of interest.
- * @param cosmo The cosmological model.
- */
-__attribute__((always_inline)) INLINE static float
-hydro_get_physical_internal_energy_dt(const struct part* restrict p,
-                                      const struct cosmology* cosmo) {
-  error("Needs implementing");
-  return 0.f;
-}
-
-/**
- * @brief Check if the gradient matrix for this particle is well behaved.
- */
-#define hydro_part_geometry_well_behaved(p) \
-  (p->geometry.wcorr > const_gizmo_min_wcorr)
-
-/**
- * @brief Macro used to access the name of the density field in the part struct.
- */
-#define hydro_part_get_density_variable() rho
-
-/**
- * @brief Macro used to access the name of the pressure field in the part
- * struct.
- */
-#define hydro_part_get_pressure_variable() P
-
 #endif /* SWIFT_GIZMO_MFM_HYDRO_GETTERS_H */
diff --git a/src/hydro/Gizmo/MFM/hydro_part.h b/src/hydro/Gizmo/MFM/hydro_part.h
index 117fd470b5e723f8db1044aca2c7ee8a107b7200..e708ee0b56fdfc1bb53ea8f2f9f5601fbb4077b6 100644
--- a/src/hydro/Gizmo/MFM/hydro_part.h
+++ b/src/hydro/Gizmo/MFM/hydro_part.h
@@ -19,39 +19,6 @@
 #ifndef SWIFT_GIZMO_MFM_HYDRO_PART_H
 #define SWIFT_GIZMO_MFM_HYDRO_PART_H
 
-#include "black_holes_struct.h"
-#include "chemistry_struct.h"
-#include "cooling_struct.h"
-#include "star_formation_struct.h"
-#include "timestep_limiter_struct.h"
-#include "tracers_struct.h"
-
-/* Extra particle data not needed during the computation. */
-struct xpart {
-
-  /* Offset between current position and position at last tree rebuild. */
-  float x_diff[3];
-
-  /* Offset between the current position and position at the last sort. */
-  float x_diff_sort[3];
-
-  /* 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;
-
-  /* Additional data used by the tracers */
-  struct tracers_xpart_data tracers_data;
-
-  /* Additional data used by the star formation */
-  struct star_formation_xpart_data sf_data;
-
-} SWIFT_STRUCT_ALIGN;
-
 /* Data of a single particle. */
 struct part {
 
@@ -64,8 +31,16 @@ struct part {
   /* Particle position. */
   double x[3];
 
-  /* Particle predicted velocity. */
-  float v[3];
+  /* In MFM, the particle and fluid velocities are the same.
+     We use an anonymous union to make sure we can reference the
+     same array with both names. */
+  union {
+    /* Particle predicted velocity. */
+    float v[3];
+
+    /* Fluid velocity. */
+    float fluid_v[3];
+  };
 
   /* Particle acceleration. */
   float a_hydro[3];
diff --git a/src/hydro/Gizmo/MFM/hydro_setters.h b/src/hydro/Gizmo/MFM/hydro_setters.h
index a503833daccdce24124dfa25c64575230e3e5327..a367f73f3dc7f76acc1f4c98ca3dac035f9d958a 100644
--- a/src/hydro/Gizmo/MFM/hydro_setters.h
+++ b/src/hydro/Gizmo/MFM/hydro_setters.h
@@ -19,54 +19,6 @@
 #ifndef SWIFT_GIZMO_MFM_HYDRO_SETTERS_H
 #define SWIFT_GIZMO_MFM_HYDRO_SETTERS_H
 
-#include "const.h"
-
-/**
- * @brief Set the primitive variables for the given particle to the given
- * values.
- *
- * @param p Particle.
- * @param W Primitive variables.
- */
-__attribute__((always_inline)) INLINE static void
-hydro_part_set_primitive_variables(struct part* restrict p, const float* W) {
-
-  p->rho = W[0];
-  p->v[0] = W[1];
-  p->v[1] = W[2];
-  p->v[2] = W[3];
-  p->P = W[4];
-}
-
-/**
- * @brief Set the conserved variables for the given particle to the given
- * values.
- *
- * @param p Particle.
- * @param Q Conserved variables.
- */
-__attribute__((always_inline)) INLINE static void
-hydro_part_set_conserved_variables(struct part* restrict p, const float* Q) {
-
-  p->conserved.mass = Q[0];
-  p->conserved.momentum[0] = Q[1];
-  p->conserved.momentum[2] = Q[2];
-  p->conserved.momentum[3] = Q[3];
-  p->conserved.energy = Q[4];
-}
-
-/**
- * @brief Set the correction value for degenerate particle configurations for
- * the given particle to the given value.
- *
- * @param p Particle.
- * @param wcorr New value.
- */
-__attribute__((always_inline)) INLINE static void hydro_part_set_wcorr(
-    struct part* restrict p, const float wcorr) {
-  p->geometry.wcorr = wcorr;
-}
-
 /**
  * @brief Reset the fluxes for the given particle.
  *
@@ -128,323 +80,4 @@ hydro_part_update_fluxes_right(struct part* restrict p, const float* fluxes,
 #endif
 }
 
-/**
- * @brief Set the gradients for the given particle to zero.
- *
- * @param p Particle.
- */
-__attribute__((always_inline)) INLINE static void hydro_part_reset_gradients(
-    struct part* restrict p) {
-
-  p->gradients.rho[0] = 0.0f;
-  p->gradients.rho[1] = 0.0f;
-  p->gradients.rho[2] = 0.0f;
-
-  p->gradients.v[0][0] = 0.0f;
-  p->gradients.v[0][1] = 0.0f;
-  p->gradients.v[0][2] = 0.0f;
-  p->gradients.v[1][0] = 0.0f;
-  p->gradients.v[1][1] = 0.0f;
-  p->gradients.v[1][2] = 0.0f;
-  p->gradients.v[2][0] = 0.0f;
-  p->gradients.v[2][1] = 0.0f;
-  p->gradients.v[2][2] = 0.0f;
-
-  p->gradients.P[0] = 0.0f;
-  p->gradients.P[1] = 0.0f;
-  p->gradients.P[2] = 0.0f;
-}
-
-/**
- * @brief Set the gradients for the given particle to the given values.
- *
- * @param p Particle.
- */
-__attribute__((always_inline)) INLINE static void hydro_part_set_gradients(
-    struct part* restrict p, const float* gradrho, const float* gradvx,
-    const float* gradvy, const float* gradvz, const float* gradP) {
-
-  p->gradients.rho[0] = gradrho[0];
-  p->gradients.rho[1] = gradrho[1];
-  p->gradients.rho[2] = gradrho[2];
-
-  p->gradients.v[0][0] = gradvx[0];
-  p->gradients.v[0][1] = gradvx[1];
-  p->gradients.v[0][2] = gradvx[2];
-  p->gradients.v[1][0] = gradvy[0];
-  p->gradients.v[1][1] = gradvy[1];
-  p->gradients.v[1][2] = gradvy[2];
-  p->gradients.v[2][0] = gradvz[0];
-  p->gradients.v[2][1] = gradvz[1];
-  p->gradients.v[2][2] = gradvz[2];
-
-  p->gradients.P[0] = gradP[0];
-  p->gradients.P[1] = gradP[1];
-  p->gradients.P[2] = gradP[2];
-}
-
-/**
- * @brief Update the gradients for the given particle with the given
- * contributions.
- *
- * @param p Particle.
- * @param drho Density gradient contribution.
- * @param dvx x velocity gradient contribution.
- * @param dvy y velocity gradient contribution.
- * @param dvz z velocity gradient contribution.
- * @param dP Pressure gradient contribution.
- */
-__attribute__((always_inline)) INLINE static void hydro_part_update_gradients(
-    struct part* restrict p, const float* drho, const float* dvx,
-    const float* dvy, const float* dvz, const float* dP) {
-
-  p->gradients.rho[0] += drho[0];
-  p->gradients.rho[1] += drho[1];
-  p->gradients.rho[2] += drho[2];
-
-  p->gradients.v[0][0] += dvx[0];
-  p->gradients.v[0][1] += dvx[1];
-  p->gradients.v[0][2] += dvx[2];
-  p->gradients.v[1][0] += dvy[0];
-  p->gradients.v[1][1] += dvy[1];
-  p->gradients.v[1][2] += dvy[2];
-  p->gradients.v[2][0] += dvz[0];
-  p->gradients.v[2][1] += dvz[1];
-  p->gradients.v[2][2] += dvz[2];
-
-  p->gradients.P[0] += dP[0];
-  p->gradients.P[1] += dP[1];
-  p->gradients.P[2] += dP[2];
-}
-
-/**
- * @brief Normalise the gradients for the given particle with the given
- * normalisation factor.
- *
- * @param p Particle.
- * @param norm Normalisation factor.
- */
-__attribute__((always_inline)) INLINE static void
-hydro_part_normalise_gradients(struct part* restrict p, const float norm) {
-
-  p->gradients.rho[0] *= norm;
-  p->gradients.rho[1] *= norm;
-  p->gradients.rho[2] *= norm;
-
-  p->gradients.v[0][0] *= norm;
-  p->gradients.v[0][1] *= norm;
-  p->gradients.v[0][2] *= norm;
-  p->gradients.v[1][0] *= norm;
-  p->gradients.v[1][1] *= norm;
-  p->gradients.v[1][2] *= norm;
-  p->gradients.v[2][0] *= norm;
-  p->gradients.v[2][1] *= norm;
-  p->gradients.v[2][2] *= norm;
-
-  p->gradients.P[0] *= norm;
-  p->gradients.P[1] *= norm;
-  p->gradients.P[2] *= norm;
-}
-
-/**
- * @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 Sets the time derivative of the co-moving internal energy of a
- * particle
- *
- * We assume a constant density for the conversion to entropy.
- *
- * @param p The particle of interest.
- * @param du_dt The new time derivative of the comoving internal energy.
- */
-__attribute__((always_inline)) INLINE static void
-hydro_set_comoving_internal_energy_dt(struct part* restrict p,
-                                      const float du_dt) {
-  error("Needs implementing");
-}
-
-/**
- * @brief Sets the time derivative of the physical internal energy of a particle
- *
- * We assume a constant density for the conversion to entropy.
- *
- * @param p The particle of interest.
- * @param cosmo Cosmology data structure
- * @param du_dt The time derivative of the physical internal energy.
- */
-__attribute__((always_inline)) INLINE static void
-hydro_set_physical_internal_energy_dt(struct part* restrict p,
-                                      const struct cosmology* restrict cosmo,
-                                      const float du_dt) {
-  error("Needs implementing");
-}
-/**
- * @brief Sets the physical entropy of a particle
- *
- * @param p The particle of interest.
- * @param xp The extended particle data.
- * @param cosmo Cosmology data structure
- * @param entropy The physical entropy
- */
-__attribute__((always_inline)) INLINE static void hydro_set_physical_entropy(
-    struct part* p, struct xpart* xp, const struct cosmology* cosmo,
-    const float entropy) {
-
-  error("Needs implementing");
-}
-
-/**
- * @brief Sets the physical internal energy of a particle
- *
- * @param p The particle of interest.
- * @param xp The extended particle data.
- * @param cosmo Cosmology data structure
- * @param u The physical internal energy
- */
-__attribute__((always_inline)) INLINE static void
-hydro_set_physical_internal_energy(struct part* p, struct xpart* xp,
-                                   const struct cosmology* cosmo,
-                                   const float u) {
-  error("Need implementing");
-}
-
-/**
- * @brief Sets the drifted physical internal energy of a particle
- *
- * @param p The particle of interest.
- * @param cosmo Cosmology data structure
- * @param u The physical internal energy
- */
-__attribute__((always_inline)) INLINE static void
-hydro_set_drifted_physical_internal_energy(struct part* p,
-                                           const struct cosmology* cosmo,
-                                           const float u) {
-  error("Need implementing");
-}
-
-/**
- * @brief Update the value of the viscosity alpha for the scheme.
- *
- * @param p the particle of interest
- * @param alpha the new value for the viscosity coefficient.
- */
-__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
-    struct part* restrict p, float alpha) {
-  /* Purposefully left empty */
-}
-
-/**
- * @brief Update the value of the viscosity alpha to the
- *        feedback reset value for the scheme.
- *
- * @param p the particle of interest
- */
-__attribute__((always_inline)) INLINE static void
-hydro_diffusive_feedback_reset(struct part* restrict p) {
-  /* Purposefully left empty */
-}
-
-/**
- * @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->rho;
-}
-
-/**
- * @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->rho;
-}
-
-/**
- * @brief Modifies the thermal state of a particle to the imposed internal
- * energy
- *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
- *
- * @param p The particle
- * @param u The new internal energy
- */
-__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
-    struct part* restrict p, float u) {
-
-  /* conserved.energy is NOT the specific energy (u), but the total thermal
-     energy (u*m) */
-  p->conserved.energy = u * p->conserved.mass;
-#ifdef GIZMO_TOTAL_ENERGY
-  /* add the kinetic energy */
-  p->conserved.energy +=
-      0.5f * p->conserved.mass *
-      (p->conserved.momentum[0] * p->v[0] + p->conserved.momentum[1] * p->v[1] +
-       p->conserved.momentum[2] * p->v[2]);
-#endif
-  p->P = hydro_gamma_minus_one * p->rho * u;
-}
-
-/**
- * @brief Modifies the thermal state of a particle to the imposed entropy
- *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
- *
- * @param p The particle
- * @param S The new entropy
- */
-__attribute__((always_inline)) INLINE static void hydro_set_entropy(
-    struct part* restrict p, float S) {
-
-  p->conserved.energy = S * pow_gamma_minus_one(p->rho) *
-                        hydro_one_over_gamma_minus_one * p->conserved.mass;
-#ifdef GIZMO_TOTAL_ENERGY
-  /* add the kinetic energy */
-  p->conserved.energy +=
-      0.5f * p->conserved.mass *
-      (p->conserved.momentum[0] * p->v[0] + p->conserved.momentum[1] * p->v[1] +
-       p->conserved.momentum[2] * p->v[2]);
-#endif
-  p->P = S * pow_gamma(p->rho);
-}
-
-/**
- * @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) {
-
-  /* We store the initial energy per unit mass in the energy
-   * variable as the conversion to energy will be done later,
-   * in hydro_first_init_part(). */
-  p->conserved.energy = u_init;
-}
-
 #endif /* SWIFT_GIZMO_MFM_HYDRO_SETTERS_H */
diff --git a/src/hydro/Gizmo/MFM/hydro_slope_limiters_cell.h b/src/hydro/Gizmo/MFM/hydro_slope_limiters_cell.h
deleted file mode 100644
index 87ceda43e02085174d998266dcbbf391d118712a..0000000000000000000000000000000000000000
--- a/src/hydro/Gizmo/MFM/hydro_slope_limiters_cell.h
+++ /dev/null
@@ -1,75 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
- *
- * 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_GIZMO_MFM_SLOPE_LIMITER_CELL_H
-#define SWIFT_GIZMO_MFM_SLOPE_LIMITER_CELL_H
-
-#include <float.h>
-
-/**
- * @brief Initialize variables for the cell wide slope limiter
- *
- * @param p Particle.
- */
-__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
-    struct part* p) {
-
-  p->limiter.rho[0] = FLT_MAX;
-  p->limiter.rho[1] = -FLT_MAX;
-  p->limiter.v[0][0] = FLT_MAX;
-  p->limiter.v[0][1] = -FLT_MAX;
-  p->limiter.v[1][0] = FLT_MAX;
-  p->limiter.v[1][1] = -FLT_MAX;
-  p->limiter.v[2][0] = FLT_MAX;
-  p->limiter.v[2][1] = -FLT_MAX;
-  p->limiter.P[0] = FLT_MAX;
-  p->limiter.P[1] = -FLT_MAX;
-
-  p->limiter.maxr = -FLT_MAX;
-}
-
-/**
- * @brief Collect information for the cell wide slope limiter during the
- * neighbour loop
- *
- * @param pi Particle i.
- * @param pj Particle j.
- * @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) {
-
-  /* basic slope limiter: collect the maximal and the minimal value for the
-   * primitive variables among the ngbs */
-  pi->limiter.rho[0] = min(pj->rho, pi->limiter.rho[0]);
-  pi->limiter.rho[1] = max(pj->rho, pi->limiter.rho[1]);
-
-  pi->limiter.v[0][0] = min(pj->v[0], pi->limiter.v[0][0]);
-  pi->limiter.v[0][1] = max(pj->v[0], pi->limiter.v[0][1]);
-  pi->limiter.v[1][0] = min(pj->v[1], pi->limiter.v[1][0]);
-  pi->limiter.v[1][1] = max(pj->v[1], pi->limiter.v[1][1]);
-  pi->limiter.v[2][0] = min(pj->v[2], pi->limiter.v[2][0]);
-  pi->limiter.v[2][1] = max(pj->v[2], pi->limiter.v[2][1]);
-
-  pi->limiter.P[0] = min(pj->P, pi->limiter.P[0]);
-  pi->limiter.P[1] = max(pj->P, pi->limiter.P[1]);
-
-  pi->limiter.maxr = max(r, pi->limiter.maxr);
-}
-
-#endif /* SWIFT_GIZMO_MFM_SLOPE_LIMITER_CELL_H */
diff --git a/src/hydro/Gizmo/MFM/hydro_velocities.h b/src/hydro/Gizmo/MFM/hydro_velocities.h
new file mode 100644
index 0000000000000000000000000000000000000000..b729197a4881e379601518c446fa4f753ceeb1f2
--- /dev/null
+++ b/src/hydro/Gizmo/MFM/hydro_velocities.h
@@ -0,0 +1,104 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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_GIZMO_MFM_HYDRO_VELOCITIES_H
+#define SWIFT_GIZMO_MFM_HYDRO_VELOCITIES_H
+
+#ifdef GIZMO_FIX_PARTICLES
+#error "Fixed particles are not allowed for GIZMO MFM!"
+#endif
+
+#ifdef GIZMO_STEER_MOTION
+#error "Steering particle movement is not allowed for GIZMO MFM!"
+#endif
+
+/**
+ * @brief Initialize the GIZMO particle velocities before the start of the
+ * actual run based on the initial value of the primitive velocity.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_velocities_init(
+    struct part* restrict p, struct xpart* restrict xp) {
+
+  xp->v_full[0] = p->v[0];
+  xp->v_full[1] = p->v[1];
+  xp->v_full[2] = p->v[2];
+}
+
+/**
+ * @brief Set the particle velocity field that will be used to deboost fluid
+ * velocities during the force loop.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_velocities_prepare_force(struct part* restrict p,
+                               const struct xpart* restrict xp) {}
+
+/**
+ * @brief Set the variables that will be used to update the smoothing length
+ * during the drift (these will depend on the movement of the particles).
+ *
+ * @param p The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_velocities_end_force(
+    struct part* restrict p) {
+
+  /* Add normalization to h_dt. */
+  p->force.h_dt *= p->h * hydro_dimension_inv;
+}
+
+/**
+ * @brief Set the velocity of a GIZMO particle, based on the values of its
+ * primitive variables and the geometry of its mesh-free "cell".
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_velocities_set(
+    struct part* restrict p, struct xpart* restrict xp) {
+
+  /* Set the velocities: */
+  /* We first set the particle velocity */
+  if (p->conserved.mass > 0.0f && p->rho > 0.0f) {
+
+    const float inverse_mass = 1.0f / p->conserved.mass;
+
+    /* Normal case: set particle velocity to fluid velocity. */
+    xp->v_full[0] = p->conserved.momentum[0] * inverse_mass;
+    xp->v_full[1] = p->conserved.momentum[1] * inverse_mass;
+    xp->v_full[2] = p->conserved.momentum[2] * inverse_mass;
+
+  } else {
+    /* Vacuum particles have no fluid velocity. */
+    xp->v_full[0] = 0.0f;
+    xp->v_full[1] = 0.0f;
+    xp->v_full[2] = 0.0f;
+  }
+
+  if (p->gpart) {
+    p->gpart->v_full[0] = xp->v_full[0];
+    p->gpart->v_full[1] = xp->v_full[1];
+    p->gpart->v_full[2] = xp->v_full[2];
+  }
+}
+
+#endif /* SWIFT_GIZMO_MFM_HYDRO_VELOCITIES_H */
diff --git a/src/hydro/Gizmo/MFV/hydro_debug.h b/src/hydro/Gizmo/MFV/hydro_debug.h
index 412fbf70a799b653735b288b67c7f91ae54a2f98..95a773712a487e50831a9f57f99620e46f266057 100644
--- a/src/hydro/Gizmo/MFV/hydro_debug.h
+++ b/src/hydro/Gizmo/MFV/hydro_debug.h
@@ -55,23 +55,18 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "wcount=%.3e}\n",
       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], p->h, p->time_bin, p->limiter_data.wakeup,
-      p->primitives.v[0], p->primitives.v[1], p->primitives.v[2],
-      p->primitives.rho, p->primitives.P, p->primitives.gradients.rho[0],
-      p->primitives.gradients.rho[1], p->primitives.gradients.rho[2],
-      p->primitives.gradients.v[0][0], p->primitives.gradients.v[0][1],
-      p->primitives.gradients.v[0][2], p->primitives.gradients.v[1][0],
-      p->primitives.gradients.v[1][1], p->primitives.gradients.v[1][2],
-      p->primitives.gradients.v[2][0], p->primitives.gradients.v[2][1],
-      p->primitives.gradients.v[2][2], p->primitives.gradients.P[0],
-      p->primitives.gradients.P[1], p->primitives.gradients.P[2],
-      p->primitives.limiter.rho[0], p->primitives.limiter.rho[1],
-      p->primitives.limiter.v[0][0], p->primitives.limiter.v[0][1],
-      p->primitives.limiter.v[1][0], p->primitives.limiter.v[1][1],
-      p->primitives.limiter.v[2][0], p->primitives.limiter.v[2][1],
-      p->primitives.limiter.P[0], p->primitives.limiter.P[1],
-      p->primitives.limiter.maxr, p->conserved.momentum[0],
-      p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.mass,
-      p->conserved.energy, p->geometry.volume, p->geometry.matrix_E[0][0],
+      p->fluid_v[0], p->fluid_v[1], p->fluid_v[2], p->rho, p->P,
+      p->gradients.rho[0], p->gradients.rho[1], p->gradients.rho[2],
+      p->gradients.v[0][0], p->gradients.v[0][1], p->gradients.v[0][2],
+      p->gradients.v[1][0], p->gradients.v[1][1], p->gradients.v[1][2],
+      p->gradients.v[2][0], p->gradients.v[2][1], p->gradients.v[2][2],
+      p->gradients.P[0], p->gradients.P[1], p->gradients.P[2],
+      p->limiter.rho[0], p->limiter.rho[1], p->limiter.v[0][0],
+      p->limiter.v[0][1], p->limiter.v[1][0], p->limiter.v[1][1],
+      p->limiter.v[2][0], p->limiter.v[2][1], p->limiter.P[0], p->limiter.P[1],
+      p->limiter.maxr, p->conserved.momentum[0], p->conserved.momentum[1],
+      p->conserved.momentum[2], p->conserved.mass, p->conserved.energy,
+      p->geometry.volume, p->geometry.matrix_E[0][0],
       p->geometry.matrix_E[0][1], p->geometry.matrix_E[0][2],
       p->geometry.matrix_E[1][0], p->geometry.matrix_E[1][1],
       p->geometry.matrix_E[1][2], p->geometry.matrix_E[2][0],
diff --git a/src/hydro/Gizmo/MFV/hydro_getters.h b/src/hydro/Gizmo/MFV/hydro_getters.h
index 811ff0e36f11b1f59c54f290d77a2ff5ff4e9653..9c8293328fd14cc4fba2d78104be3fb9b8eabb77 100644
--- a/src/hydro/Gizmo/MFV/hydro_getters.h
+++ b/src/hydro/Gizmo/MFV/hydro_getters.h
@@ -19,120 +19,6 @@
 #ifndef SWIFT_GIZMO_MFV_HYDRO_GETTERS_H
 #define SWIFT_GIZMO_MFV_HYDRO_GETTERS_H
 
-#include "cosmology.h"
-#include "equation_of_state.h"
-
-/**
- * @brief Get a 5-element state vector W containing the primitive hydrodynamic
- * variables.
- *
- * @param p Particle.
- * @param W Pointer to the array in which the result needs to be stored (of size
- * 5 or more).
- */
-__attribute__((always_inline)) INLINE static void
-hydro_part_get_primitive_variables(const struct part* restrict p, float* W) {
-
-  W[0] = p->primitives.rho;
-  W[1] = p->primitives.v[0];
-  W[2] = p->primitives.v[1];
-  W[3] = p->primitives.v[2];
-  W[4] = p->primitives.P;
-}
-
-/**
- * @brief Get the degenerate case correction factor for the given particle.
- *
- * @param p Particle.
- * @return Degenerate case correction factor.
- */
-__attribute__((always_inline)) INLINE static float hydro_part_get_wcorr(
-    const struct part* restrict p) {
-  return p->density.wcorr;
-}
-
-/**
- * @brief Get the gradients of the primitive variables for the given particle.
- *
- * @param p Particle.
- * @param drho Density gradient (of size 3 or more).
- * @param ddvx x velocity gradient (of size 3 or more).
- * @param ddvy y velocity gradient (of size 3 or more).
- * @param ddvz z velocity gradient (of size 3 or more).
- * @param dP Pressure gradient (of size 3 or more).
- */
-__attribute__((always_inline)) INLINE static void hydro_part_get_gradients(
-    const struct part* restrict p, float* drho, float* dvx, float* dvy,
-    float* dvz, float* dP) {
-
-  drho[0] = p->primitives.gradients.rho[0];
-  drho[1] = p->primitives.gradients.rho[1];
-  drho[2] = p->primitives.gradients.rho[2];
-
-  dvx[0] = p->primitives.gradients.v[0][0];
-  dvx[1] = p->primitives.gradients.v[0][1];
-  dvx[2] = p->primitives.gradients.v[0][2];
-  dvy[0] = p->primitives.gradients.v[1][0];
-  dvy[1] = p->primitives.gradients.v[1][1];
-  dvy[2] = p->primitives.gradients.v[1][2];
-  dvz[0] = p->primitives.gradients.v[2][0];
-  dvz[1] = p->primitives.gradients.v[2][1];
-  dvz[2] = p->primitives.gradients.v[2][2];
-
-  dP[0] = p->primitives.gradients.P[0];
-  dP[1] = p->primitives.gradients.P[1];
-  dP[2] = p->primitives.gradients.P[2];
-}
-
-/**
- * @brief Get the pressure gradient for the given particle.
- *
- * @param p Particle.
- * @param gradP Pressure gradient (of size 3 or more).
- */
-__attribute__((always_inline)) INLINE static void
-hydro_part_get_pressure_gradient(const struct part* restrict p,
-                                 float gradP[3]) {
-
-  gradP[0] = p->primitives.gradients.P[0];
-  gradP[1] = p->primitives.gradients.P[1];
-  gradP[2] = p->primitives.gradients.P[2];
-}
-
-/**
- * @brief Get the slope limiter variables for the given particle.
- *
- * @param p Particle.
- * @param rholim Minimum and maximum density of neighbours (of size 2 or more).
- * @param vxlim Minimum and maximum x velocity of neighbours (of size 2 or
- * more).
- * @param vylim Minimum and maximum y velocity of neighbours (of size 2 or
- * more).
- * @param vzlim Minimum and maximum z velocity of neighbours (of size 2 or
- * more).
- * @param Plim Minimum and maximum pressure of neighbours (of size 2 or more).
- * @param rmax Maximum distance of any neighbour (of size 1 or more).
- */
-__attribute__((always_inline)) INLINE static void hydro_part_get_slope_limiter(
-    const struct part* restrict p, float* rholim, float* vxlim, float* vylim,
-    float* vzlim, float* Plim, float* rmax) {
-
-  rholim[0] = p->primitives.limiter.rho[0];
-  rholim[1] = p->primitives.limiter.rho[1];
-
-  vxlim[0] = p->primitives.limiter.v[0][0];
-  vxlim[1] = p->primitives.limiter.v[0][1];
-  vylim[0] = p->primitives.limiter.v[1][0];
-  vylim[1] = p->primitives.limiter.v[1][1];
-  vzlim[0] = p->primitives.limiter.v[2][0];
-  vzlim[1] = p->primitives.limiter.v[2][1];
-
-  Plim[0] = p->primitives.limiter.P[0];
-  Plim[1] = p->primitives.limiter.P[1];
-
-  rmax[0] = p->primitives.limiter.maxr;
-}
-
 /**
  * @brief Get the fluxes for the given particle.
  *
@@ -142,240 +28,11 @@ __attribute__((always_inline)) INLINE static void hydro_part_get_slope_limiter(
 __attribute__((always_inline)) INLINE static void hydro_part_get_fluxes(
     const struct part* restrict p, float* flux) {
 
-  flux[0] = p->conserved.flux.mass;
-  flux[1] = p->conserved.flux.momentum[0];
-  flux[2] = p->conserved.flux.momentum[1];
-  flux[3] = p->conserved.flux.momentum[2];
-  flux[4] = p->conserved.flux.energy;
-}
-
-/**
- * @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 physical internal energy of a particle
- *
- * @param p The particle of interest.
- * @param xp The extended data of 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 xpart* restrict xp,
-                                   const struct cosmology* cosmo) {
-
-  return cosmo->a_factor_internal_energy *
-         hydro_get_comoving_internal_energy(p);
+  flux[0] = p->flux.mass;
+  flux[1] = p->flux.momentum[0];
+  flux[2] = p->flux.momentum[1];
+  flux[3] = p->flux.momentum[2];
+  flux[4] = p->flux.energy;
 }
 
-/**
- * @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_drifted_physical_internal_energy(const struct part* restrict p,
-                                           const struct cosmology* cosmo) {
-
-  return hydro_get_physical_internal_energy(p, /*xp=*/NULL, cosmo);
-}
-
-/**
- * @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 physical internal energy of a particle
- *
- * @param p The particle of interest.
- * @param xp The extended data of 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 xpart* restrict xp,
-    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 internal energy of a particle
- *
- * @param p The particle of interest.
- * @param cosmo The cosmological model.
- */
-__attribute__((always_inline)) INLINE static float
-hydro_get_drifted_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 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 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
- */
-__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
-    const struct part* restrict p) {
-
-  return p->primitives.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 mass of a particle
- *
- * @param p The particle of interest
- */
-__attribute__((always_inline)) INLINE static float hydro_get_mass(
-    const struct part* restrict p) {
-
-  return p->conserved.mass;
-}
-
-/**
- * @brief Returns the velocities drifted to the current time of a particle.
- *
- * @param p The particle of interest
- * @param xp The extended data of the particle.
- * @param dt_kick_hydro The time (for hydro accelerations) since the last kick.
- * @param dt_kick_grav The time (for gravity accelerations) since the last kick.
- * @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_kick_hydro,
-    float dt_kick_grav, float v[3]) {
-
-  if (p->conserved.mass > 0.) {
-    v[0] = p->primitives.v[0] +
-           p->conserved.flux.momentum[0] * dt_kick_hydro / p->conserved.mass;
-    v[1] = p->primitives.v[1] +
-           p->conserved.flux.momentum[1] * dt_kick_hydro / p->conserved.mass;
-    v[2] = p->primitives.v[2] +
-           p->conserved.flux.momentum[2] * dt_kick_hydro / p->conserved.mass;
-  } else {
-    v[0] = p->primitives.v[0];
-    v[1] = p->primitives.v[1];
-    v[2] = p->primitives.v[2];
-  }
-
-  // MATTHIEU: Bert is this correct?
-  v[0] += xp->a_grav[0] * dt_kick_grav;
-  v[1] += xp->a_grav[1] * dt_kick_grav;
-  v[2] += xp->a_grav[2] * dt_kick_grav;
-}
-
-/**
- * @brief Returns the time derivative of co-moving internal energy of a particle
- *
- * We assume a constant density.
- *
- * @param p The particle of interest
- */
-__attribute__((always_inline)) INLINE static float
-hydro_get_comoving_internal_energy_dt(const struct part* restrict p) {
-
-  error("Needs implementing");
-  return 0.f;
-}
-
-/**
- * @brief Returns the time derivative of physical internal energy of a particle
- *
- * We assume a constant density.
- *
- * @param p The particle of interest.
- * @param cosmo The cosmological model.
- */
-__attribute__((always_inline)) INLINE static float
-hydro_get_physical_internal_energy_dt(const struct part* restrict p,
-                                      const struct cosmology* cosmo) {
-  error("Needs implementing");
-  return 0.f;
-}
-
-/**
- * @brief Check if the gradient matrix for this particle is well behaved.
- */
-#define hydro_part_geometry_well_behaved(p) \
-  (p->density.wcorr > const_gizmo_min_wcorr)
-
-/**
- * @brief Macro used to access the name of the density field in the part struct.
- */
-#define hydro_part_get_density_variable() primitives.rho
-
-/**
- * @brief Macro used to access the name of the pressure field in the part
- * struct.
- */
-#define hydro_part_get_pressure_variable() primitives.P
-
 #endif /* SWIFT_GIZMO_MFV_HYDRO_GETTERS_H */
diff --git a/src/hydro/Gizmo/MFV/hydro_part.h b/src/hydro/Gizmo/MFV/hydro_part.h
index fe210ea24ea6858878855e60fdc9e1959558534d..82a7296c17f21ef4c0cf5416b498201d78a965d3 100644
--- a/src/hydro/Gizmo/MFV/hydro_part.h
+++ b/src/hydro/Gizmo/MFV/hydro_part.h
@@ -19,39 +19,6 @@
 #ifndef SWIFT_GIZMO_MFV_HYDRO_PART_H
 #define SWIFT_GIZMO_MFV_HYDRO_PART_H
 
-#include "black_holes_struct.h"
-#include "chemistry_struct.h"
-#include "cooling_struct.h"
-#include "star_formation_struct.h"
-#include "timestep_limiter_struct.h"
-#include "tracers_struct.h"
-
-/* Extra particle data not needed during the computation. */
-struct xpart {
-
-  /* Offset between current position and position at last tree rebuild. */
-  float x_diff[3];
-
-  /* Offset between the current position and position at the last sort. */
-  float x_diff_sort[3];
-
-  /* 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;
-
-  /* Additional data used by the tracers */
-  struct tracers_xpart_data tracers_data;
-
-  /* Additional data used by the star formation */
-  struct star_formation_xpart_data sf_data;
-
-} SWIFT_STRUCT_ALIGN;
-
 /* Data of a single particle. */
 struct part {
 
@@ -73,50 +40,45 @@ struct part {
   /* Particle smoothing length. */
   float h;
 
-  /* The primitive hydrodynamical variables. */
-  struct {
-
-    /* Density. */
-    float rho;
-
-    /* Fluid velocity. */
-    float v[3];
+  /* Density. */
+  float rho;
 
-    /* Pressure. */
-    float P;
+  /* Fluid velocity. */
+  float fluid_v[3];
 
-    /* Gradients of the primitive variables. */
-    struct {
+  /* Pressure. */
+  float P;
 
-      /* Density gradients. */
-      float rho[3];
+  /* Gradients of the primitive variables. */
+  struct {
 
-      /* Fluid velocity gradients. */
-      float v[3][3];
+    /* Density gradients. */
+    float rho[3];
 
-      /* Pressure gradients. */
-      float P[3];
+    /* Fluid velocity gradients. */
+    float v[3][3];
 
-    } gradients;
+    /* Pressure gradients. */
+    float P[3];
 
-    /* Quantities needed by the slope limiter. */
-    struct {
+  } gradients;
 
-      /* Extreme values of the density among the neighbours. */
-      float rho[2];
+  /* Quantities needed by the slope limiter. */
+  struct {
 
-      /* Extreme values of the fluid velocity among the neighbours. */
-      float v[3][2];
+    /* Extreme values of the density among the neighbours. */
+    float rho[2];
 
-      /* Extreme values of the pressure among the neighbours. */
-      float P[2];
+    /* Extreme values of the fluid velocity among the neighbours. */
+    float v[3][2];
 
-      /* Maximal distance to all neighbouring faces. */
-      float maxr;
+    /* Extreme values of the pressure among the neighbours. */
+    float P[2];
 
-    } limiter;
+    /* Maximal distance to all neighbouring faces. */
+    float maxr;
 
-  } primitives;
+  } limiter;
 
   /* The conserved hydrodynamical variables. */
   struct {
@@ -130,21 +92,21 @@ struct part {
     /* Fluid thermal energy (not per unit mass!). */
     float energy;
 
-    /* Fluxes. */
-    struct {
+  } conserved;
 
-      /* Mass flux. */
-      float mass;
+  /* Fluxes. */
+  struct {
 
-      /* Momentum flux. */
-      float momentum[3];
+    /* Mass flux. */
+    float mass;
 
-      /* Energy flux. */
-      float energy;
+    /* Momentum flux. */
+    float momentum[3];
 
-    } flux;
+    /* Energy flux. */
+    float energy;
 
-  } conserved;
+  } flux;
 
   /* Geometrical quantities used for hydro. */
   struct {
@@ -159,6 +121,9 @@ struct part {
     /* Centroid of the "cell". */
     float centroid[3];
 
+    /* Correction factor for wcount. */
+    float wcorr;
+
   } geometry;
 
   /* Variables used for timestep calculation. */
@@ -181,9 +146,6 @@ struct part {
     /* Particle number density. */
     float wcount;
 
-    /* Correction factor for wcount. */
-    float wcorr;
-
   } density;
 
   /* Quantities used during the force loop. */
diff --git a/src/hydro/Gizmo/MFV/hydro_setters.h b/src/hydro/Gizmo/MFV/hydro_setters.h
index 770336ea63ee325b71fdbd5c2648bc9e2ff9189c..f37a1eebc7b870773d7880b4f631a5eb931f3577 100644
--- a/src/hydro/Gizmo/MFV/hydro_setters.h
+++ b/src/hydro/Gizmo/MFV/hydro_setters.h
@@ -19,54 +19,6 @@
 #ifndef SWIFT_GIZMO_MFV_HYDRO_SETTERS_H
 #define SWIFT_GIZMO_MFV_HYDRO_SETTERS_H
 
-#include "const.h"
-
-/**
- * @brief Set the primitive variables for the given particle to the given
- * values.
- *
- * @param p Particle.
- * @param W Primitive variables.
- */
-__attribute__((always_inline)) INLINE static void
-hydro_part_set_primitive_variables(struct part* restrict p, const float* W) {
-
-  p->primitives.rho = W[0];
-  p->primitives.v[0] = W[1];
-  p->primitives.v[1] = W[2];
-  p->primitives.v[2] = W[3];
-  p->primitives.P = W[4];
-}
-
-/**
- * @brief Set the conserved variables for the given particle to the given
- * values.
- *
- * @param p Particle.
- * @param Q Conserved variables.
- */
-__attribute__((always_inline)) INLINE static void
-hydro_part_set_conserved_variables(struct part* restrict p, const float* Q) {
-
-  p->conserved.mass = Q[0];
-  p->conserved.momentum[0] = Q[1];
-  p->conserved.momentum[2] = Q[2];
-  p->conserved.momentum[3] = Q[3];
-  p->conserved.energy = Q[4];
-}
-
-/**
- * @brief Set the correction value for degenerate particle configurations for
- * the given particle to the given value.
- *
- * @param p Particle.
- * @param wcorr New value.
- */
-__attribute__((always_inline)) INLINE static void hydro_part_set_wcorr(
-    struct part* restrict p, const float wcorr) {
-  p->density.wcorr = wcorr;
-}
-
 /**
  * @brief Reset the fluxes for the given particle.
  *
@@ -75,11 +27,11 @@ __attribute__((always_inline)) INLINE static void hydro_part_set_wcorr(
 __attribute__((always_inline)) INLINE static void hydro_part_reset_fluxes(
     struct part* restrict p) {
 
-  p->conserved.flux.mass = 0.0f;
-  p->conserved.flux.momentum[0] = 0.0f;
-  p->conserved.flux.momentum[1] = 0.0f;
-  p->conserved.flux.momentum[2] = 0.0f;
-  p->conserved.flux.energy = 0.0f;
+  p->flux.mass = 0.0f;
+  p->flux.momentum[0] = 0.0f;
+  p->flux.momentum[1] = 0.0f;
+  p->flux.momentum[2] = 0.0f;
+  p->flux.energy = 0.0f;
 
   p->gravity.mflux[0] = 0.0f;
   p->gravity.mflux[1] = 0.0f;
@@ -101,20 +53,20 @@ __attribute__((always_inline)) INLINE static void hydro_part_update_fluxes_left(
   p->gravity.mflux[1] += fluxes[0] * dx[1];
   p->gravity.mflux[2] += fluxes[0] * dx[2];
 
-  p->conserved.flux.mass -= fluxes[0];
-  p->conserved.flux.momentum[0] -= fluxes[1];
-  p->conserved.flux.momentum[1] -= fluxes[2];
-  p->conserved.flux.momentum[2] -= fluxes[3];
-  p->conserved.flux.energy -= fluxes[4];
+  p->flux.mass -= fluxes[0];
+  p->flux.momentum[0] -= fluxes[1];
+  p->flux.momentum[1] -= fluxes[2];
+  p->flux.momentum[2] -= fluxes[3];
+  p->flux.energy -= fluxes[4];
 
 #ifndef GIZMO_TOTAL_ENERGY
-  const float ekin = 0.5f * (p->primitives.v[0] * p->primitives.v[0] +
-                             p->primitives.v[1] * p->primitives.v[1] +
-                             p->primitives.v[2] * p->primitives.v[2]);
-  p->conserved.flux.energy += fluxes[1] * p->primitives.v[0];
-  p->conserved.flux.energy += fluxes[2] * p->primitives.v[1];
-  p->conserved.flux.energy += fluxes[3] * p->primitives.v[2];
-  p->conserved.flux.energy -= fluxes[0] * ekin;
+  const float ekin =
+      0.5f * (p->fluid_v[0] * p->fluid_v[0] + p->fluid_v[1] * p->fluid_v[1] +
+              p->fluid_v[2] * p->fluid_v[2]);
+  p->flux.energy += fluxes[1] * p->fluid_v[0];
+  p->flux.energy += fluxes[2] * p->fluid_v[1];
+  p->flux.energy += fluxes[3] * p->fluid_v[2];
+  p->flux.energy -= fluxes[0] * ekin;
 #endif
 }
 
@@ -134,340 +86,21 @@ hydro_part_update_fluxes_right(struct part* restrict p, const float* fluxes,
   p->gravity.mflux[1] += fluxes[0] * dx[1];
   p->gravity.mflux[2] += fluxes[0] * dx[2];
 
-  p->conserved.flux.mass += fluxes[0];
-  p->conserved.flux.momentum[0] += fluxes[1];
-  p->conserved.flux.momentum[1] += fluxes[2];
-  p->conserved.flux.momentum[2] += fluxes[3];
-  p->conserved.flux.energy += fluxes[4];
+  p->flux.mass += fluxes[0];
+  p->flux.momentum[0] += fluxes[1];
+  p->flux.momentum[1] += fluxes[2];
+  p->flux.momentum[2] += fluxes[3];
+  p->flux.energy += fluxes[4];
 
 #ifndef GIZMO_TOTAL_ENERGY
-  const float ekin = 0.5f * (p->primitives.v[0] * p->primitives.v[0] +
-                             p->primitives.v[1] * p->primitives.v[1] +
-                             p->primitives.v[2] * p->primitives.v[2]);
-  p->conserved.flux.energy -= fluxes[1] * p->primitives.v[0];
-  p->conserved.flux.energy -= fluxes[2] * p->primitives.v[1];
-  p->conserved.flux.energy -= fluxes[3] * p->primitives.v[2];
-  p->conserved.flux.energy += fluxes[0] * ekin;
+  const float ekin =
+      0.5f * (p->fluid_v[0] * p->fluid_v[0] + p->fluid_v[1] * p->fluid_v[1] +
+              p->fluid_v[2] * p->fluid_v[2]);
+  p->flux.energy -= fluxes[1] * p->fluid_v[0];
+  p->flux.energy -= fluxes[2] * p->fluid_v[1];
+  p->flux.energy -= fluxes[3] * p->fluid_v[2];
+  p->flux.energy += fluxes[0] * ekin;
 #endif
 }
 
-/**
- * @brief Set the gradients for the given particle to zero.
- *
- * @param p Particle.
- */
-__attribute__((always_inline)) INLINE static void hydro_part_reset_gradients(
-    struct part* restrict p) {
-
-  p->primitives.gradients.rho[0] = 0.0f;
-  p->primitives.gradients.rho[1] = 0.0f;
-  p->primitives.gradients.rho[2] = 0.0f;
-
-  p->primitives.gradients.v[0][0] = 0.0f;
-  p->primitives.gradients.v[0][1] = 0.0f;
-  p->primitives.gradients.v[0][2] = 0.0f;
-  p->primitives.gradients.v[1][0] = 0.0f;
-  p->primitives.gradients.v[1][1] = 0.0f;
-  p->primitives.gradients.v[1][2] = 0.0f;
-  p->primitives.gradients.v[2][0] = 0.0f;
-  p->primitives.gradients.v[2][1] = 0.0f;
-  p->primitives.gradients.v[2][2] = 0.0f;
-
-  p->primitives.gradients.P[0] = 0.0f;
-  p->primitives.gradients.P[1] = 0.0f;
-  p->primitives.gradients.P[2] = 0.0f;
-}
-
-/**
- * @brief Set the gradients for the given particle to the given values.
- *
- * @param p Particle.
- */
-__attribute__((always_inline)) INLINE static void hydro_part_set_gradients(
-    struct part* restrict p, const float* gradrho, const float* gradvx,
-    const float* gradvy, const float* gradvz, const float* gradP) {
-
-  p->primitives.gradients.rho[0] = gradrho[0];
-  p->primitives.gradients.rho[1] = gradrho[1];
-  p->primitives.gradients.rho[2] = gradrho[2];
-
-  p->primitives.gradients.v[0][0] = gradvx[0];
-  p->primitives.gradients.v[0][1] = gradvx[1];
-  p->primitives.gradients.v[0][2] = gradvx[2];
-  p->primitives.gradients.v[1][0] = gradvy[0];
-  p->primitives.gradients.v[1][1] = gradvy[1];
-  p->primitives.gradients.v[1][2] = gradvy[2];
-  p->primitives.gradients.v[2][0] = gradvz[0];
-  p->primitives.gradients.v[2][1] = gradvz[1];
-  p->primitives.gradients.v[2][2] = gradvz[2];
-
-  p->primitives.gradients.P[0] = gradP[0];
-  p->primitives.gradients.P[1] = gradP[1];
-  p->primitives.gradients.P[2] = gradP[2];
-}
-
-/**
- * @brief Update the gradients for the given particle with the given
- * contributions.
- *
- * @param p Particle.
- * @param drho Density gradient contribution.
- * @param dvx x velocity gradient contribution.
- * @param dvy y velocity gradient contribution.
- * @param dvz z velocity gradient contribution.
- * @param dP Pressure gradient contribution.
- */
-__attribute__((always_inline)) INLINE static void hydro_part_update_gradients(
-    struct part* restrict p, const float* drho, const float* dvx,
-    const float* dvy, const float* dvz, const float* dP) {
-
-  p->primitives.gradients.rho[0] += drho[0];
-  p->primitives.gradients.rho[1] += drho[1];
-  p->primitives.gradients.rho[2] += drho[2];
-
-  p->primitives.gradients.v[0][0] += dvx[0];
-  p->primitives.gradients.v[0][1] += dvx[1];
-  p->primitives.gradients.v[0][2] += dvx[2];
-  p->primitives.gradients.v[1][0] += dvy[0];
-  p->primitives.gradients.v[1][1] += dvy[1];
-  p->primitives.gradients.v[1][2] += dvy[2];
-  p->primitives.gradients.v[2][0] += dvz[0];
-  p->primitives.gradients.v[2][1] += dvz[1];
-  p->primitives.gradients.v[2][2] += dvz[2];
-
-  p->primitives.gradients.P[0] += dP[0];
-  p->primitives.gradients.P[1] += dP[1];
-  p->primitives.gradients.P[2] += dP[2];
-}
-
-/**
- * @brief Normalise the gradients for the given particle with the given
- * normalisation factor.
- *
- * @param p Particle.
- * @param norm Normalisation factor.
- */
-__attribute__((always_inline)) INLINE static void
-hydro_part_normalise_gradients(struct part* restrict p, const float norm) {
-
-  p->primitives.gradients.rho[0] *= norm;
-  p->primitives.gradients.rho[1] *= norm;
-  p->primitives.gradients.rho[2] *= norm;
-
-  p->primitives.gradients.v[0][0] *= norm;
-  p->primitives.gradients.v[0][1] *= norm;
-  p->primitives.gradients.v[0][2] *= norm;
-  p->primitives.gradients.v[1][0] *= norm;
-  p->primitives.gradients.v[1][1] *= norm;
-  p->primitives.gradients.v[1][2] *= norm;
-  p->primitives.gradients.v[2][0] *= norm;
-  p->primitives.gradients.v[2][1] *= norm;
-  p->primitives.gradients.v[2][2] *= norm;
-
-  p->primitives.gradients.P[0] *= norm;
-  p->primitives.gradients.P[1] *= norm;
-  p->primitives.gradients.P[2] *= norm;
-}
-
-/**
- * @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 Sets the time derivative of the co-moving internal energy of a
- * particle
- *
- * We assume a constant density for the conversion to entropy.
- *
- * @param p The particle of interest.
- * @param du_dt The new time derivative of the comoving internal energy.
- */
-__attribute__((always_inline)) INLINE static void
-hydro_set_comoving_internal_energy_dt(struct part* restrict p,
-                                      const float du_dt) {
-  error("Needs implementing");
-}
-
-/**
- * @brief Sets the time derivative of the physical internal energy of a particle
- *
- * We assume a constant density for the conversion to entropy.
- *
- * @param p The particle of interest.
- * @param cosmo Cosmology data structure
- * @param du_dt The time derivative of the physical internal energy.
- */
-__attribute__((always_inline)) INLINE static void
-hydro_set_physical_internal_energy_dt(struct part* restrict p,
-                                      const struct cosmology* restrict cosmo,
-                                      const float du_dt) {
-  error("Needs implementing");
-}
-/**
- * @brief Sets the physical entropy of a particle
- *
- * @param p The particle of interest.
- * @param xp The extended particle data.
- * @param cosmo Cosmology data structure
- * @param entropy The physical entropy
- */
-__attribute__((always_inline)) INLINE static void hydro_set_physical_entropy(
-    struct part* p, struct xpart* xp, const struct cosmology* cosmo,
-    const float entropy) {
-
-  error("Needs implementing");
-}
-
-/**
- * @brief Sets the physical internal energy of a particle
- *
- * @param p The particle of interest.
- * @param xp The extended particle data.
- * @param cosmo Cosmology data structure
- * @param u The physical internal energy
- */
-__attribute__((always_inline)) INLINE static void
-hydro_set_physical_internal_energy(struct part* p, struct xpart* xp,
-                                   const struct cosmology* cosmo,
-                                   const float u) {
-  error("Need implementing");
-}
-
-/**
- * @brief Sets the drifted physical internal energy of a particle
- *
- * @param p The particle of interest.
- * @param cosmo Cosmology data structure
- * @param u The physical internal energy
- */
-__attribute__((always_inline)) INLINE static void
-hydro_set_drifted_physical_internal_energy(struct part* p,
-                                           const struct cosmology* cosmo,
-                                           const float u) {
-  error("Need implementing");
-}
-
-/**
- * @brief Update the value of the viscosity alpha for the scheme.
- *
- * @param p the particle of interest
- * @param alpha the new value for the viscosity coefficient.
- */
-__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
-    struct part* restrict p, float alpha) {
-  /* Purposefully left empty */
-}
-
-/**
- * @brief Update the value of the viscosity alpha to the
- *        feedback reset value for the scheme.
- *
- * @param p the particle of interest
- */
-__attribute__((always_inline)) INLINE static void
-hydro_diffusive_feedback_reset(struct part* restrict p) {
-  /* Purposefully left empty */
-}
-
-/**
- * @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 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;
-}
-
-/**
- * @brief Modifies the thermal state of a particle to the imposed internal
- * energy
- *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
- *
- * @param p The particle
- * @param u The new internal energy
- */
-__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
-    struct part* restrict p, float u) {
-
-  /* conserved.energy is NOT the specific energy (u), but the total thermal
-     energy (u*m) */
-  p->conserved.energy = u * 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;
-}
-
-/**
- * @brief Modifies the thermal state of a particle to the imposed entropy
- *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
- *
- * @param p The particle
- * @param S The new entropy
- */
-__attribute__((always_inline)) INLINE static void hydro_set_entropy(
-    struct part* restrict p, float S) {
-
-  p->conserved.energy = S * pow_gamma_minus_one(p->primitives.rho) *
-                        hydro_one_over_gamma_minus_one * 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 = S * pow_gamma(p->primitives.rho);
-}
-
-/**
- * @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) {
-
-  /* We store the initial energy per unit mass in the energy
-   * variable as the conversion to energy will be done later,
-   * in hydro_first_init_part(). */
-  p->conserved.energy = u_init;
-}
-
 #endif /* SWIFT_GIZMO_MFV_HYDRO_SETTERS_H */
diff --git a/src/hydro/Gizmo/MFV/hydro_slope_limiters_cell.h b/src/hydro/Gizmo/MFV/hydro_slope_limiters_cell.h
deleted file mode 100644
index d270ff15f31a0ce124eaabaf9aaf9c8dd469cd9a..0000000000000000000000000000000000000000
--- a/src/hydro/Gizmo/MFV/hydro_slope_limiters_cell.h
+++ /dev/null
@@ -1,85 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
- *
- * 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_GIZMO_MFV_SLOPE_LIMITER_CELL_H
-#define SWIFT_GIZMO_MFV_SLOPE_LIMITER_CELL_H
-
-#include <float.h>
-
-/**
- * @brief Initialize variables for the cell wide slope limiter
- *
- * @param p Particle.
- */
-__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
-    struct part* p) {
-
-  p->primitives.limiter.rho[0] = FLT_MAX;
-  p->primitives.limiter.rho[1] = -FLT_MAX;
-  p->primitives.limiter.v[0][0] = FLT_MAX;
-  p->primitives.limiter.v[0][1] = -FLT_MAX;
-  p->primitives.limiter.v[1][0] = FLT_MAX;
-  p->primitives.limiter.v[1][1] = -FLT_MAX;
-  p->primitives.limiter.v[2][0] = FLT_MAX;
-  p->primitives.limiter.v[2][1] = -FLT_MAX;
-  p->primitives.limiter.P[0] = FLT_MAX;
-  p->primitives.limiter.P[1] = -FLT_MAX;
-
-  p->primitives.limiter.maxr = -FLT_MAX;
-}
-
-/**
- * @brief Collect information for the cell wide slope limiter during the
- * neighbour loop
- *
- * @param pi Particle i.
- * @param pj Particle j.
- * @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) {
-
-  /* basic slope limiter: collect the maximal and the minimal value for the
-   * primitive variables among the ngbs */
-  pi->primitives.limiter.rho[0] =
-      min(pj->primitives.rho, pi->primitives.limiter.rho[0]);
-  pi->primitives.limiter.rho[1] =
-      max(pj->primitives.rho, pi->primitives.limiter.rho[1]);
-
-  pi->primitives.limiter.v[0][0] =
-      min(pj->primitives.v[0], pi->primitives.limiter.v[0][0]);
-  pi->primitives.limiter.v[0][1] =
-      max(pj->primitives.v[0], pi->primitives.limiter.v[0][1]);
-  pi->primitives.limiter.v[1][0] =
-      min(pj->primitives.v[1], pi->primitives.limiter.v[1][0]);
-  pi->primitives.limiter.v[1][1] =
-      max(pj->primitives.v[1], pi->primitives.limiter.v[1][1]);
-  pi->primitives.limiter.v[2][0] =
-      min(pj->primitives.v[2], pi->primitives.limiter.v[2][0]);
-  pi->primitives.limiter.v[2][1] =
-      max(pj->primitives.v[2], pi->primitives.limiter.v[2][1]);
-
-  pi->primitives.limiter.P[0] =
-      min(pj->primitives.P, pi->primitives.limiter.P[0]);
-  pi->primitives.limiter.P[1] =
-      max(pj->primitives.P, pi->primitives.limiter.P[1]);
-
-  pi->primitives.limiter.maxr = max(r, pi->primitives.limiter.maxr);
-}
-
-#endif /* SWIFT_GIZMO_MFV_SLOPE_LIMITER_CELL_H */
diff --git a/src/hydro/Gizmo/MFV/hydro_velocities.h b/src/hydro/Gizmo/MFV/hydro_velocities.h
index a61a482683291cfc0662116ada99b94a4ccfe68d..189c27b19bf98396e939098044554c2075723d29 100644
--- a/src/hydro/Gizmo/MFV/hydro_velocities.h
+++ b/src/hydro/Gizmo/MFV/hydro_velocities.h
@@ -16,8 +16,8 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef SWIFT_HYDRO_VELOCITIES_H
-#define SWIFT_HYDRO_VELOCITIES_H
+#ifndef SWIFT_GIZMO_MFV_HYDRO_VELOCITIES_H
+#define SWIFT_GIZMO_MFV_HYDRO_VELOCITIES_H
 
 /**
  * @brief Initialize the GIZMO particle velocities before the start of the
@@ -30,13 +30,13 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_init(
     struct part* restrict p, struct xpart* restrict xp) {
 
 #ifdef GIZMO_FIX_PARTICLES
-  p->v[0] = 0.f;
-  p->v[1] = 0.f;
-  p->v[2] = 0.f;
+  p->v[0] = 0.0f;
+  p->v[1] = 0.0f;
+  p->v[2] = 0.0f;
 #else
-  p->v[0] = p->primitives.v[0];
-  p->v[1] = p->primitives.v[1];
-  p->v[2] = p->primitives.v[2];
+  p->v[0] = p->fluid_v[0];
+  p->v[1] = p->fluid_v[1];
+  p->v[2] = p->fluid_v[2];
 #endif
 
   xp->v_full[0] = p->v[0];
@@ -87,15 +87,15 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_set(
 /* We first set the particle velocity. */
 #ifdef GIZMO_FIX_PARTICLES
 
-  p->v[0] = 0.f;
-  p->v[1] = 0.f;
-  p->v[2] = 0.f;
+  p->v[0] = 0.0f;
+  p->v[1] = 0.0f;
+  p->v[2] = 0.0f;
 
 #else  // GIZMO_FIX_PARTICLES
 
-  if (p->conserved.mass > 0.f && p->primitives.rho > 0.f) {
+  if (p->conserved.mass > 0.0f && p->rho > 0.0f) {
 
-    const float inverse_mass = 1.f / p->conserved.mass;
+    const float inverse_mass = 1.0f / p->conserved.mass;
 
     /* Normal case: set particle velocity to fluid velocity. */
     p->v[0] = p->conserved.momentum[0] * inverse_mass;
@@ -116,15 +116,14 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_set(
     const float R = get_radius_dimension_sphere(p->geometry.volume);
     const float eta = 0.25f;
     const float etaR = eta * R;
-    const float xi = 1.f;
-    const float soundspeed =
-        sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+    const float xi = 1.0f;
+    const float soundspeed = sqrtf(hydro_gamma * p->P / p->rho);
     /* We only apply the correction if the offset between centroid and position
        is too large. */
     if (d > 0.9f * etaR) {
       float fac = xi * soundspeed / d;
       if (d < 1.1f * etaR) {
-        fac *= 5.f * (d - 0.9f * etaR) / etaR;
+        fac *= 5.0f * (d - 0.9f * etaR) / etaR;
       }
       p->v[0] -= ds[0] * fac;
       p->v[1] -= ds[1] * fac;
@@ -134,9 +133,9 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_set(
 #endif  // GIZMO_STEER_MOTION
   } else {
     /* Vacuum particles have no fluid velocity. */
-    p->v[0] = 0.f;
-    p->v[1] = 0.f;
-    p->v[2] = 0.f;
+    p->v[0] = 0.0f;
+    p->v[1] = 0.0f;
+    p->v[2] = 0.0f;
   }
 
 #endif  // GIZMO_FIX_PARTICLES
@@ -153,4 +152,4 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_set(
   }
 }
 
-#endif /* SWIFT_HYDRO_VELOCITIES_H */
+#endif /* SWIFT_GIZMO_MFV_HYDRO_VELOCITIES_H */
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 23bd906e78942cb6f58b0a67a2d55c5a8cb2062d..ba4ce84ef1673b7de44c30eb60019d584185c113 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -27,10 +27,7 @@
 #include "hydro_gradients.h"
 #include "hydro_setters.h"
 #include "hydro_space.h"
-
-#if defined(GIZMO_MFV_SPH)
-#include "MFV/hydro_velocities.h"
-#endif
+#include "hydro_velocities.h"
 
 #include <float.h>
 
@@ -172,15 +169,8 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
   hydro_part_set_primitive_variables(p, W);
   hydro_part_set_conserved_variables(p, Q);
 
-#if defined(GIZMO_MFV_SPH)
   /* initialize the particle velocity based on the primitive fluid velocity */
   hydro_velocities_init(p, xp);
-#elif defined(GIZMO_MFM_SPH)
-  /* initialize the particle velocity based on the primitive fluid velocity */
-  xp->v_full[0] = p->v[0];
-  xp->v_full[1] = p->v[1];
-  xp->v_full[2] = p->v[2];
-#endif
 
   /* ignore accelerations present in the initial condition */
   p->a_hydro[0] = 0.0f;
@@ -451,10 +441,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
 
   hydro_gradients_init(p);
 
-#if defined(GIZMO_MFV_SPH)
-  /* Set the actual velocity of the particle */
   hydro_velocities_prepare_force(p, xp);
-#endif
 }
 
 /**
@@ -669,13 +656,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 __attribute__((always_inline)) INLINE static void hydro_end_force(
     struct part* p, const struct cosmology* cosmo) {
 
-  /* Add normalization to h_dt. */
-  p->force.h_dt *= p->h * hydro_dimension_inv;
-
-  // MATTHIEU: Bert is this correct? Do we need cosmology terms here?
-#if defined(GIZMO_MFV_SPH)
   hydro_velocities_end_force(p);
-#endif
 }
 
 /**
@@ -769,11 +750,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
       hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
   if (p->conserved.energy < min_energy * p->conserved.mass) {
     p->conserved.energy = min_energy * p->conserved.mass;
-#if defined(GIZMO_MFV_SPH)
-    p->conserved.flux.energy = 0.0f;
-#elif defined(GIZMO_MFM_SPH)
     p->flux.energy = 0.0f;
-#endif
   }
 
   // MATTHIEU: Apply the entropy floor here.
@@ -804,33 +781,9 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     /* Make sure the gpart knows the mass has changed. */
     p->gpart->mass = p->conserved.mass;
   }
+#endif
 
   hydro_velocities_set(p, xp);
-#elif defined(GIZMO_MFM_SPH)
-  /* Set the velocities: */
-  /* We first set the particle velocity */
-  if (p->conserved.mass > 0.0f && p->rho > 0.0f) {
-
-    const float inverse_mass = 1.0f / p->conserved.mass;
-
-    /* Normal case: set particle velocity to fluid velocity. */
-    xp->v_full[0] = p->conserved.momentum[0] * inverse_mass;
-    xp->v_full[1] = p->conserved.momentum[1] * inverse_mass;
-    xp->v_full[2] = p->conserved.momentum[2] * inverse_mass;
-
-  } else {
-    /* Vacuum particles have no fluid velocity. */
-    xp->v_full[0] = 0.0f;
-    xp->v_full[1] = 0.0f;
-    xp->v_full[2] = 0.0f;
-  }
-
-  if (p->gpart) {
-    p->gpart->v_full[0] = xp->v_full[0];
-    p->gpart->v_full[1] = xp->v_full[1];
-    p->gpart->v_full[2] = xp->v_full[2];
-  }
-#endif
 
 #ifdef GIZMO_LLOYD_ITERATION
   /* reset conserved variables to safe values */
diff --git a/src/hydro/Gizmo/hydro_getters.h b/src/hydro/Gizmo/hydro_getters.h
index 59217744441c00680b20ed9ab13f35ba7a22c5cb..a7ed402c1fbefdaf9e4acaf4c2b6db53a7b9b4e3 100644
--- a/src/hydro/Gizmo/hydro_getters.h
+++ b/src/hydro/Gizmo/hydro_getters.h
@@ -19,10 +19,345 @@
 #ifndef SWIFT_GIZMO_HYDRO_GETTERS_H
 #define SWIFT_GIZMO_HYDRO_GETTERS_H
 
+#include "cosmology.h"
+#include "equation_of_state.h"
+
 #if defined(GIZMO_MFV_SPH)
 #include "MFV/hydro_getters.h"
 #elif defined(GIZMO_MFM_SPH)
 #include "MFM/hydro_getters.h"
 #endif
 
+/**
+ * @brief Get a 5-element state vector W containing the primitive hydrodynamic
+ * variables.
+ *
+ * @param p Particle.
+ * @param W Pointer to the array in which the result needs to be stored (of size
+ * 5 or more).
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_part_get_primitive_variables(const struct part* restrict p, float* W) {
+
+  W[0] = p->rho;
+  W[1] = p->fluid_v[0];
+  W[2] = p->fluid_v[1];
+  W[3] = p->fluid_v[2];
+  W[4] = p->P;
+}
+
+/**
+ * @brief Get the degenerate case correction factor for the given particle.
+ *
+ * @param p Particle.
+ * @return Degenerate case correction factor.
+ */
+__attribute__((always_inline)) INLINE static float hydro_part_get_wcorr(
+    const struct part* restrict p) {
+  return p->geometry.wcorr;
+}
+
+/**
+ * @brief Get the gradients of the primitive variables for the given particle.
+ *
+ * @param p Particle.
+ * @param drho Density gradient (of size 3 or more).
+ * @param ddvx x velocity gradient (of size 3 or more).
+ * @param ddvy y velocity gradient (of size 3 or more).
+ * @param ddvz z velocity gradient (of size 3 or more).
+ * @param dP Pressure gradient (of size 3 or more).
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_get_gradients(
+    const struct part* restrict p, float* drho, float* dvx, float* dvy,
+    float* dvz, float* dP) {
+
+  drho[0] = p->gradients.rho[0];
+  drho[1] = p->gradients.rho[1];
+  drho[2] = p->gradients.rho[2];
+
+  dvx[0] = p->gradients.v[0][0];
+  dvx[1] = p->gradients.v[0][1];
+  dvx[2] = p->gradients.v[0][2];
+  dvy[0] = p->gradients.v[1][0];
+  dvy[1] = p->gradients.v[1][1];
+  dvy[2] = p->gradients.v[1][2];
+  dvz[0] = p->gradients.v[2][0];
+  dvz[1] = p->gradients.v[2][1];
+  dvz[2] = p->gradients.v[2][2];
+
+  dP[0] = p->gradients.P[0];
+  dP[1] = p->gradients.P[1];
+  dP[2] = p->gradients.P[2];
+}
+
+/**
+ * @brief Get the pressure gradient for the given particle.
+ *
+ * @param p Particle.
+ * @param gradP Pressure gradient (of size 3 or more).
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_part_get_pressure_gradient(const struct part* restrict p,
+                                 float gradP[3]) {
+
+  gradP[0] = p->gradients.P[0];
+  gradP[1] = p->gradients.P[1];
+  gradP[2] = p->gradients.P[2];
+}
+
+/**
+ * @brief Get the slope limiter variables for the given particle.
+ *
+ * @param p Particle.
+ * @param rholim Minimum and maximum density of neighbours (of size 2 or more).
+ * @param vxlim Minimum and maximum x velocity of neighbours (of size 2 or
+ * more).
+ * @param vylim Minimum and maximum y velocity of neighbours (of size 2 or
+ * more).
+ * @param vzlim Minimum and maximum z velocity of neighbours (of size 2 or
+ * more).
+ * @param Plim Minimum and maximum pressure of neighbours (of size 2 or more).
+ * @param rmax Maximum distance of any neighbour (of size 1 or more).
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_get_slope_limiter(
+    const struct part* restrict p, float* rholim, float* vxlim, float* vylim,
+    float* vzlim, float* Plim, float* rmax) {
+
+  rholim[0] = p->limiter.rho[0];
+  rholim[1] = p->limiter.rho[1];
+
+  vxlim[0] = p->limiter.v[0][0];
+  vxlim[1] = p->limiter.v[0][1];
+  vylim[0] = p->limiter.v[1][0];
+  vylim[1] = p->limiter.v[1][1];
+  vzlim[0] = p->limiter.v[2][0];
+  vzlim[1] = p->limiter.v[2][1];
+
+  Plim[0] = p->limiter.P[0];
+  Plim[1] = p->limiter.P[1];
+
+  rmax[0] = p->limiter.maxr;
+}
+
+/**
+ * @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->rho > 0.0f)
+    return gas_internal_energy_from_pressure(p->rho, p->P);
+  else
+    return 0.;
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param xp The extended data of 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 xpart* restrict xp,
+                                   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_drifted_physical_internal_energy(const struct part* restrict p,
+                                           const struct cosmology* cosmo) {
+
+  return hydro_get_physical_internal_energy(p, /*xp=*/NULL, cosmo);
+}
+
+/**
+ * @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->rho > 0.0f) {
+    return gas_entropy_from_pressure(p->rho, p->P);
+  } else {
+    return 0.;
+  }
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param xp The extended data of 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 xpart* restrict xp,
+    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 internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_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 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->rho > 0.0f)
+    return gas_soundspeed_from_pressure(p->rho, p->P);
+  else
+    return 0.;
+}
+
+/**
+ * @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
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
+    const struct part* restrict p) {
+
+  return p->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->P;
+}
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_mass(
+    const struct part* restrict p) {
+
+  return p->conserved.mass;
+}
+
+/**
+ * @brief Returns the velocities drifted to the current time of a particle.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle.
+ * @param dt_kick_hydro The time (for hydro accelerations) since the last kick.
+ * @param dt_kick_grav The time (for gravity accelerations) since the last kick.
+ * @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_kick_hydro,
+    float dt_kick_grav, float v[3]) {
+
+  if (p->conserved.mass > 0.) {
+    const float m_inv = 1.0f / p->conserved.mass;
+    v[0] = p->fluid_v[0] + p->flux.momentum[0] * dt_kick_hydro * m_inv;
+    v[1] = p->fluid_v[1] + p->flux.momentum[1] * dt_kick_hydro * m_inv;
+    v[2] = p->fluid_v[2] + p->flux.momentum[2] * dt_kick_hydro * m_inv;
+  } else {
+    v[0] = p->fluid_v[0];
+    v[1] = p->fluid_v[1];
+    v[2] = p->fluid_v[2];
+  }
+
+  // MATTHIEU: Bert is this correct?
+  v[0] += xp->a_grav[0] * dt_kick_grav;
+  v[1] += xp->a_grav[1] * dt_kick_grav;
+  v[2] += xp->a_grav[2] * dt_kick_grav;
+}
+
+/**
+ * @brief Returns the time derivative of co-moving internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy_dt(const struct part* restrict p) {
+
+  error("Needs implementing");
+  return 0.0f;
+}
+
+/**
+ * @brief Returns the time derivative of physical internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy_dt(const struct part* restrict p,
+                                      const struct cosmology* cosmo) {
+  error("Needs implementing");
+  return 0.0f;
+}
+
+/**
+ * @brief Check if the gradient matrix for this particle is well behaved.
+ *
+ * @param p Particle.
+ * @return 1 if the gradient matrix is well behaved, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int
+hydro_part_geometry_well_behaved(const struct part* restrict p) {
+
+  return p->geometry.wcorr > const_gizmo_min_wcorr;
+}
+
 #endif /* SWIFT_GIZMO_HYDRO_GETTERS_H */
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index bb15093600cfe67fd74109d6e0ac038875cae05a..3fc8d78ef0bddb67262173ea0f17d2ee962efb40 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -21,7 +21,6 @@
 
 #include "adiabatic_index.h"
 #include "hydro.h"
-#include "hydro_getters.h"
 #include "hydro_gradients.h"
 #include "hydro_slope_limiters.h"
 #include "io_properties.h"
@@ -63,9 +62,8 @@ INLINE static void hydro_read_particles(struct part* parts,
                                 UNIT_CONV_NO_UNITS, parts, id);
   list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
                                 UNIT_CONV_ACCELERATION, parts, a_hydro);
-  list[7] =
-      io_make_input_field("Density", FLOAT, 1, OPTIONAL, UNIT_CONV_DENSITY,
-                          parts, hydro_part_get_density_variable());
+  list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
+                                UNIT_CONV_DENSITY, parts, rho);
 }
 
 /**
@@ -218,7 +216,7 @@ INLINE static void hydro_write_particles(const struct part* parts,
                            parts, id, "Unique IDs of the particles");
 
   list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
-                                 parts, hydro_part_get_density_variable(),
+                                 parts, rho,
                                  "Co-moving mass densities of the particles");
 
   list[7] = io_make_output_field_convert_part(
@@ -226,8 +224,7 @@ INLINE static void hydro_write_particles(const struct part* parts,
       "Co-moving entropies of the particles");
 
   list[8] = io_make_output_field("Pressures", FLOAT, 1, UNIT_CONV_PRESSURE,
-                                 -3.f * hydro_gamma, parts,
-                                 hydro_part_get_pressure_variable(),
+                                 -3.f * hydro_gamma, parts, P,
                                  "Co-moving pressures of the particles");
 
   list[9] = io_make_output_field_convert_part(
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index 8df75767909bb8d6a4a2b8ba80a9d65cd7674adf..7003396da33c80b2001433a80c287125ab0dbc58 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/Gizmo/hydro_part.h
@@ -19,6 +19,39 @@
 #ifndef SWIFT_GIZMO_HYDRO_PART_H
 #define SWIFT_GIZMO_HYDRO_PART_H
 
+#include "black_holes_struct.h"
+#include "chemistry_struct.h"
+#include "cooling_struct.h"
+#include "star_formation_struct.h"
+#include "timestep_limiter_struct.h"
+#include "tracers_struct.h"
+
+/* Extra particle data not needed during the computation. */
+struct xpart {
+
+  /* Offset between current position and position at last tree rebuild. */
+  float x_diff[3];
+
+  /* Offset between the current position and position at the last sort. */
+  float x_diff_sort[3];
+
+  /* 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;
+
+  /* Additional data used by the tracers */
+  struct tracers_xpart_data tracers_data;
+
+  /* Additional data used by the star formation */
+  struct star_formation_xpart_data sf_data;
+
+} SWIFT_STRUCT_ALIGN;
+
 /* Import the right hydro particle definition */
 #if defined(GIZMO_MFV_SPH)
 #include "MFV/hydro_part.h"
diff --git a/src/hydro/Gizmo/hydro_setters.h b/src/hydro/Gizmo/hydro_setters.h
index 351796979f84417647719dd2076b9f32fc1f551b..8fa345fce50c77d861e57be67be7bbe6fe7ad91e 100644
--- a/src/hydro/Gizmo/hydro_setters.h
+++ b/src/hydro/Gizmo/hydro_setters.h
@@ -25,4 +25,369 @@
 #include "MFM/hydro_setters.h"
 #endif
 
+/**
+ * @brief Set the primitive variables for the given particle to the given
+ * values.
+ *
+ * @param p Particle.
+ * @param W Primitive variables.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_part_set_primitive_variables(struct part* restrict p, const float* W) {
+
+  p->rho = W[0];
+  p->fluid_v[0] = W[1];
+  p->fluid_v[1] = W[2];
+  p->fluid_v[2] = W[3];
+  p->P = W[4];
+}
+
+/**
+ * @brief Set the conserved variables for the given particle to the given
+ * values.
+ *
+ * @param p Particle.
+ * @param Q Conserved variables.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_part_set_conserved_variables(struct part* restrict p, const float* Q) {
+
+  p->conserved.mass = Q[0];
+  p->conserved.momentum[0] = Q[1];
+  p->conserved.momentum[2] = Q[2];
+  p->conserved.momentum[3] = Q[3];
+  p->conserved.energy = Q[4];
+}
+
+/**
+ * @brief Set the correction value for degenerate particle configurations for
+ * the given particle to the given value.
+ *
+ * @param p Particle.
+ * @param wcorr New value.
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_set_wcorr(
+    struct part* restrict p, const float wcorr) {
+  p->geometry.wcorr = wcorr;
+}
+
+/**
+ * @brief Set the gradients for the given particle to zero.
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_reset_gradients(
+    struct part* restrict p) {
+
+  p->gradients.rho[0] = 0.0f;
+  p->gradients.rho[1] = 0.0f;
+  p->gradients.rho[2] = 0.0f;
+
+  p->gradients.v[0][0] = 0.0f;
+  p->gradients.v[0][1] = 0.0f;
+  p->gradients.v[0][2] = 0.0f;
+  p->gradients.v[1][0] = 0.0f;
+  p->gradients.v[1][1] = 0.0f;
+  p->gradients.v[1][2] = 0.0f;
+  p->gradients.v[2][0] = 0.0f;
+  p->gradients.v[2][1] = 0.0f;
+  p->gradients.v[2][2] = 0.0f;
+
+  p->gradients.P[0] = 0.0f;
+  p->gradients.P[1] = 0.0f;
+  p->gradients.P[2] = 0.0f;
+}
+
+/**
+ * @brief Set the gradients for the given particle to the given values.
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_set_gradients(
+    struct part* restrict p, const float* gradrho, const float* gradvx,
+    const float* gradvy, const float* gradvz, const float* gradP) {
+
+  p->gradients.rho[0] = gradrho[0];
+  p->gradients.rho[1] = gradrho[1];
+  p->gradients.rho[2] = gradrho[2];
+
+  p->gradients.v[0][0] = gradvx[0];
+  p->gradients.v[0][1] = gradvx[1];
+  p->gradients.v[0][2] = gradvx[2];
+  p->gradients.v[1][0] = gradvy[0];
+  p->gradients.v[1][1] = gradvy[1];
+  p->gradients.v[1][2] = gradvy[2];
+  p->gradients.v[2][0] = gradvz[0];
+  p->gradients.v[2][1] = gradvz[1];
+  p->gradients.v[2][2] = gradvz[2];
+
+  p->gradients.P[0] = gradP[0];
+  p->gradients.P[1] = gradP[1];
+  p->gradients.P[2] = gradP[2];
+}
+
+/**
+ * @brief Update the gradients for the given particle with the given
+ * contributions.
+ *
+ * @param p Particle.
+ * @param drho Density gradient contribution.
+ * @param dvx x velocity gradient contribution.
+ * @param dvy y velocity gradient contribution.
+ * @param dvz z velocity gradient contribution.
+ * @param dP Pressure gradient contribution.
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_update_gradients(
+    struct part* restrict p, const float* drho, const float* dvx,
+    const float* dvy, const float* dvz, const float* dP) {
+
+  p->gradients.rho[0] += drho[0];
+  p->gradients.rho[1] += drho[1];
+  p->gradients.rho[2] += drho[2];
+
+  p->gradients.v[0][0] += dvx[0];
+  p->gradients.v[0][1] += dvx[1];
+  p->gradients.v[0][2] += dvx[2];
+  p->gradients.v[1][0] += dvy[0];
+  p->gradients.v[1][1] += dvy[1];
+  p->gradients.v[1][2] += dvy[2];
+  p->gradients.v[2][0] += dvz[0];
+  p->gradients.v[2][1] += dvz[1];
+  p->gradients.v[2][2] += dvz[2];
+
+  p->gradients.P[0] += dP[0];
+  p->gradients.P[1] += dP[1];
+  p->gradients.P[2] += dP[2];
+}
+
+/**
+ * @brief Normalise the gradients for the given particle with the given
+ * normalisation factor.
+ *
+ * @param p Particle.
+ * @param norm Normalisation factor.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_part_normalise_gradients(struct part* restrict p, const float norm) {
+
+  p->gradients.rho[0] *= norm;
+  p->gradients.rho[1] *= norm;
+  p->gradients.rho[2] *= norm;
+
+  p->gradients.v[0][0] *= norm;
+  p->gradients.v[0][1] *= norm;
+  p->gradients.v[0][2] *= norm;
+  p->gradients.v[1][0] *= norm;
+  p->gradients.v[1][1] *= norm;
+  p->gradients.v[1][2] *= norm;
+  p->gradients.v[2][0] *= norm;
+  p->gradients.v[2][1] *= norm;
+  p->gradients.v[2][2] *= norm;
+
+  p->gradients.P[0] *= norm;
+  p->gradients.P[1] *= norm;
+  p->gradients.P[2] *= norm;
+}
+
+/**
+ * @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 Sets the time derivative of the co-moving internal energy of a
+ * particle
+ *
+ * We assume a constant density for the conversion to entropy.
+ *
+ * @param p The particle of interest.
+ * @param du_dt The new time derivative of the comoving internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_comoving_internal_energy_dt(struct part* restrict p,
+                                      const float du_dt) {
+  error("Needs implementing");
+}
+
+/**
+ * @brief Sets the time derivative of the physical internal energy of a particle
+ *
+ * We assume a constant density for the conversion to entropy.
+ *
+ * @param p The particle of interest.
+ * @param cosmo Cosmology data structure
+ * @param du_dt The time derivative of the physical internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_physical_internal_energy_dt(struct part* restrict p,
+                                      const struct cosmology* restrict cosmo,
+                                      const float du_dt) {
+  error("Needs implementing");
+}
+/**
+ * @brief Sets the physical entropy of a particle
+ *
+ * @param p The particle of interest.
+ * @param xp The extended particle data.
+ * @param cosmo Cosmology data structure
+ * @param entropy The physical entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_physical_entropy(
+    struct part* p, struct xpart* xp, const struct cosmology* cosmo,
+    const float entropy) {
+
+  error("Needs implementing");
+}
+
+/**
+ * @brief Sets the physical internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param xp The extended particle data.
+ * @param cosmo Cosmology data structure
+ * @param u The physical internal energy
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_physical_internal_energy(struct part* p, struct xpart* xp,
+                                   const struct cosmology* cosmo,
+                                   const float u) {
+  error("Need implementing");
+}
+
+/**
+ * @brief Sets the drifted physical internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param cosmo Cosmology data structure
+ * @param u The physical internal energy
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_drifted_physical_internal_energy(struct part* p,
+                                           const struct cosmology* cosmo,
+                                           const float u) {
+  error("Need implementing");
+}
+
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+    struct part* restrict p, float alpha) {
+  /* Purposefully left empty */
+}
+
+/**
+ * @brief Update the value of the viscosity alpha to the
+ *        feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_diffusive_feedback_reset(struct part* restrict p) {
+  /* Purposefully left empty */
+}
+
+/**
+ * @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->rho;
+}
+
+/**
+ * @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->rho;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed internal
+ * energy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param u The new internal energy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
+    struct part* restrict p, float u) {
+
+  /* conserved.energy is NOT the specific energy (u), but the total thermal
+     energy (u*m) */
+  p->conserved.energy = u * p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+  /* add the kinetic energy */
+  p->conserved.energy += 0.5f * p->conserved.mass *
+                         (p->conserved.momentum[0] * p->fluid_v[0] +
+                          p->conserved.momentum[1] * p->fluid_v[1] +
+                          p->conserved.momentum[2] * p->fluid_v[2]);
+#endif
+  p->P = hydro_gamma_minus_one * p->rho * u;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed entropy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param S The new entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_entropy(
+    struct part* restrict p, float S) {
+
+  p->conserved.energy = S * pow_gamma_minus_one(p->rho) *
+                        hydro_one_over_gamma_minus_one * p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+  /* add the kinetic energy */
+  p->conserved.energy += 0.5f * p->conserved.mass *
+                         (p->conserved.momentum[0] * p->fluid_v[0] +
+                          p->conserved.momentum[1] * p->fluid_v[1] +
+                          p->conserved.momentum[2] * p->fluid_v[2]);
+#endif
+  p->P = S * pow_gamma(p->rho);
+}
+
+/**
+ * @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) {
+
+  /* We store the initial energy per unit mass in the energy
+   * variable as the conversion to energy will be done later,
+   * in hydro_first_init_part(). */
+  p->conserved.energy = u_init;
+}
+
 #endif /* SWIFT_GIZMO_HYDRO_SETTERS_H */
diff --git a/src/hydro/Gizmo/hydro_slope_limiters.h b/src/hydro/Gizmo/hydro_slope_limiters.h
index bb369aa621d22c7a63312b3e73a691408a50e5bf..67f60da5c803c76a6d64c5d87dcc00b755676c6d 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters.h
+++ b/src/hydro/Gizmo/hydro_slope_limiters.h
@@ -57,12 +57,6 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
 #define HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION \
   "Cell wide slope limiter (Springel 2010)"
 
-#if defined(GIZMO_MFV_SPH)
-#include "MFV/hydro_slope_limiters_cell.h"
-#elif defined(GIZMO_MFM_SPH)
-#include "MFM/hydro_slope_limiters_cell.h"
-#endif
-
 #include "hydro_slope_limiters_cell.h"
 
 #else
diff --git a/src/hydro/Gizmo/hydro_slope_limiters_cell.h b/src/hydro/Gizmo/hydro_slope_limiters_cell.h
index b31d745dd4ba62c2959a21c5b5813d881eb7d03f..a47da9e005429e180e8813db33499463cb855425 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters_cell.h
+++ b/src/hydro/Gizmo/hydro_slope_limiters_cell.h
@@ -23,6 +23,57 @@
 #include "hydro_getters.h"
 #include "hydro_setters.h"
 
+/**
+ * @brief Initialize variables for the cell wide slope limiter
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
+    struct part* p) {
+
+  p->limiter.rho[0] = FLT_MAX;
+  p->limiter.rho[1] = -FLT_MAX;
+  p->limiter.v[0][0] = FLT_MAX;
+  p->limiter.v[0][1] = -FLT_MAX;
+  p->limiter.v[1][0] = FLT_MAX;
+  p->limiter.v[1][1] = -FLT_MAX;
+  p->limiter.v[2][0] = FLT_MAX;
+  p->limiter.v[2][1] = -FLT_MAX;
+  p->limiter.P[0] = FLT_MAX;
+  p->limiter.P[1] = -FLT_MAX;
+
+  p->limiter.maxr = -FLT_MAX;
+}
+
+/**
+ * @brief Collect information for the cell wide slope limiter during the
+ * neighbour loop
+ *
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @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) {
+
+  /* basic slope limiter: collect the maximal and the minimal value for the
+   * primitive variables among the ngbs */
+  pi->limiter.rho[0] = min(pj->rho, pi->limiter.rho[0]);
+  pi->limiter.rho[1] = max(pj->rho, pi->limiter.rho[1]);
+
+  pi->limiter.v[0][0] = min(pj->fluid_v[0], pi->limiter.v[0][0]);
+  pi->limiter.v[0][1] = max(pj->fluid_v[0], pi->limiter.v[0][1]);
+  pi->limiter.v[1][0] = min(pj->fluid_v[1], pi->limiter.v[1][0]);
+  pi->limiter.v[1][1] = max(pj->fluid_v[1], pi->limiter.v[1][1]);
+  pi->limiter.v[2][0] = min(pj->fluid_v[2], pi->limiter.v[2][0]);
+  pi->limiter.v[2][1] = max(pj->fluid_v[2], pi->limiter.v[2][1]);
+
+  pi->limiter.P[0] = min(pj->P, pi->limiter.P[0]);
+  pi->limiter.P[1] = max(pj->P, pi->limiter.P[1]);
+
+  pi->limiter.maxr = max(r, pi->limiter.maxr);
+}
+
 /**
  * @brief Slope-limit the given quantity.
  */
diff --git a/src/hydro/Gizmo/hydro_velocities.h b/src/hydro/Gizmo/hydro_velocities.h
new file mode 100644
index 0000000000000000000000000000000000000000..bb385f31728ac95ac478a1ff15852edf388a7d5e
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_velocities.h
@@ -0,0 +1,28 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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_GIZMO_HYDRO_VELOCITIES_H
+#define SWIFT_GIZMO_HYDRO_VELOCITIES_H
+
+#if defined(GIZMO_MFV_SPH)
+#include "MFV/hydro_velocities.h"
+#elif defined(GIZMO_MFM_SPH)
+#include "MFM/hydro_velocities.h"
+#endif
+
+#endif /* SWIFT_GIZMO_HYDRO_VELOCITIES_H */