diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index 2558454c60d999945a61e242f2751dd25b976d84..e00852ac2f9a04a188220075f4eb04e395f61eb6 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published
@@ -20,15 +20,17 @@
 #define SWIFT_PRESSURE_ENERGY_HYDRO_H
 
 /**
- * @file PressureEnergy/hydro.h
- * @brief Pressure-Energy implementation of SPH (Non-neighbour loop
+ * @file Minimal/hydro.h
+ * @brief Minimal conservative implementation of SPH (Non-neighbour loop
  * equations)
  *
- * The thermal variable is the energy (u) and the pressure is smoothed over
- * contact discontinuities to prevent spurious surface tension.
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
  *
- * Follows equations (16), (17) and (18) of Hopkins, P., MNRAS, 2013,
- * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
+ * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational
+ * Physics, 2012, Volume 231, Issue 3, pp. 759-794.
  */
 
 #include "adiabatic_index.h"
@@ -37,6 +39,7 @@
 #include "dimension.h"
 #include "equation_of_state.h"
 #include "hydro_properties.h"
+#include "hydro_space.h"
 #include "kernel_hydro.h"
 #include "minmax.h"
 
@@ -45,10 +48,14 @@
 /**
  * @brief Returns the comoving internal energy of a particle
  *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable.
+ *
  * @param p The particle of interest
  */
-__attribute__((always_inline)) INLINE static float hydro_get_comoving_internal_energy(
-    const struct part *restrict p) {
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy(const struct part *restrict p) {
 
   return p->u;
 }
@@ -56,11 +63,17 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_internal_e
 /**
  * @brief Returns the physical internal energy of a particle
  *
- * @param p The particle of interest
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable and converts it to
+ * physical coordinates.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
  */
-__attribute__((always_inline)) INLINE static float hydro_get_physical_internal_energy(
-    const struct part *restrict p,
-    const struct cosmology *cosmo) {
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy(const struct part *restrict p,
+                                   const struct cosmology *cosmo) {
 
   return p->u * cosmo->a_factor_internal_energy;
 }
@@ -68,6 +81,8 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_internal_e
 /**
  * @brief Returns the comoving pressure of a particle
  *
+ * Computes the pressure based on the particle's properties.
+ *
  * @param p The particle of interest
  */
 __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
@@ -79,20 +94,25 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
 /**
  * @brief Returns the physical pressure of a particle
  *
+ * Computes the pressure based on the particle's properties and
+ * convert it to physical coordinates.
+ *
  * @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) {
+    const struct part *restrict p, const struct cosmology *cosmo) {
 
-  return gas_pressure_from_internal_energy(p->rho, p->u) * cosmo->a_factor_pressure;
+  return cosmo->a_factor_pressure *
+         gas_pressure_from_internal_energy(p->rho, p->u);
 }
 
 /**
  * @brief Returns the comoving entropy of a particle
  *
- * This function actually calculates the entropy from the internal energy of the
- * particle, as in PressureEnergy (unsurprisingly) we follow U.
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable.
  *
  * @param p The particle of interest
  */
@@ -105,18 +125,19 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
 /**
  * @brief Returns the physical entropy of a particle
  *
- * This function actually calculates the entropy from the internal energy of the
- * particle, as in PressureEnergy (unsurprisingly) we follow U.
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable and converts it to
+ * physical coordinates.
  *
  * @param p The particle of interest
+ * @param cosmo The cosmological model.
  */
 __attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
-    const struct part *restrict p,
-    const struct cosmology *cosmo) {
-
-  /* Note we choose co-ordinates such that this always works out to 
-   * require no conversion */
+    const struct part *restrict p, const struct cosmology *cosmo) {
 
+  /* Note: no cosmological conversion required here with our choice of
+   * coordinates. */
   return gas_entropy_from_internal_energy(p->rho, p->u);
 }
 
@@ -125,8 +146,8 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
  *
  * @param p The particle of interest
  */
-__attribute__((always_inline)) INLINE static float hydro_get_comoving_soundspeed(
-    const struct part *restrict p) {
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_soundspeed(const struct part *restrict p) {
 
   return p->force.soundspeed;
 }
@@ -135,12 +156,13 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_soundspeed
  * @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) {
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_soundspeed(const struct part *restrict p,
+                              const struct cosmology *cosmo) {
 
-  return p->force.soundspeed * cosmo->a_factor_sound_speed;
+  return cosmo->a_factor_sound_speed * p->force.soundspeed;
 }
 
 /**
@@ -155,14 +177,15 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
 }
 
 /**
- * @brief Returns the physical density of a particle
+ * @brief Returns the comoving 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 *restrict cosmo) {
+    const struct part *restrict p, const struct cosmology *cosmo) {
 
-  return p->rho * cosmo->a3_inv;
+  return cosmo->a3_inv * p->rho;
 }
 
 /**
@@ -176,13 +199,13 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
   return p->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 The time since the last kick.
+ * @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(
@@ -190,11 +213,11 @@ __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
     float dt_kick_grav, float v[3]) {
 
   v[0] = xp->v_full[0] + p->a_hydro[0] * dt_kick_hydro +
-	                 xp->a_grav[0] * dt_kick_grav;
+         xp->a_grav[0] * dt_kick_grav;
   v[1] = xp->v_full[1] + p->a_hydro[1] * dt_kick_hydro +
-	                 xp->a_grav[1] * dt_kick_grav;
+         xp->a_grav[1] * dt_kick_grav;
   v[2] = xp->v_full[2] + p->a_hydro[2] * dt_kick_hydro +
-	                 xp->a_grav[2] * dt_kick_grav;
+         xp->a_grav[2] * dt_kick_grav;
 }
 
 /**
@@ -204,52 +227,48 @@ __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
  *
  * @param p The particle of interest
  */
-__attribute__((always_inline)) INLINE static float hydro_get_energy_dt(
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
     const struct part *restrict p) {
-  /* Since we track energy, we can forget about entropy conversions etc. */
+
   return p->u_dt;
 }
 
 /**
- * @brief Sets the time derivative of internal energy of a particle
+ * @brief Returns the time derivative of internal energy of a particle
  *
  * We assume a constant density.
  *
  * @param p The particle of interest.
  * @param du_dt The new time derivative of the internal energy.
  */
-__attribute__((always_inline)) INLINE static void hydro_set_energy_dt(
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
     struct part *restrict p, float du_dt) {
-  /* Since we track energy, we can forget about entropy conversions etc. */
+
   p->u_dt = du_dt;
 }
-
 /**
  * @brief Computes the hydro time-step of a given particle
  *
+ * This function returns the time-step of a particle given its hydro-dynamical
+ * state. A typical time-step calculation would be the use of the CFL condition.
+ *
  * @param p Pointer to the particle data
  * @param xp Pointer to the extended particle data
- *
+ * @param hydro_properties The SPH parameters
+ * @param cosmo The cosmological model.
  */
 __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
     const struct part *restrict p, const struct xpart *restrict xp,
     const struct hydro_props *restrict hydro_properties,
     const struct cosmology *restrict cosmo) {
-  
+
   const float CFL_condition = hydro_properties->CFL_condition;
 
   /* CFL condition */
-  const float dt_cfl =
-      2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
-      (p->force.v_sig * cosmo->a_factor_sound_speed);
+  const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
+                       (cosmo->a_factor_sound_speed * p->force.v_sig);
 
-  /* U change limited */
-  /* TODO: Ensure this is consistent with cosmology */
-  const float dt_u_changes = 
-    (p->u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->u_dt)
-                      : FLT_MAX;
-
-  return min(dt_u_changes, dt_cfl);
+  return dt_cfl;
 }
 
 /**
@@ -266,23 +285,21 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
  * @brief Prepares a particle for the density calculation.
  *
  * Zeroes all the relevant arrays in preparation for the sums taking place in
- * the variaous density tasks
+ * the various density loop over neighbours. Typically, all fields of the
+ * density sub-structure of a particle get zeroed in here.
  *
  * @param p The particle to act upon
+ * @param hs #hydro_space containing hydro specific space information.
  */
 __attribute__((always_inline)) INLINE static void hydro_init_part(
     struct part *restrict p, const struct hydro_space *hs) {
 
-  p->rho = 0.f;
-  p->pressure_bar = 0.f;
-  p->pressure_bar_dh = 0.f;
   p->density.wcount = 0.f;
   p->density.wcount_dh = 0.f;
-
-  p->density.div_v = 0.f;
-  p->density.rot_v[0] = 0.f;
-  p->density.rot_v[1] = 0.f;
-  p->density.rot_v[2] = 0.f;
+  p->rho = 0.f;
+  p->density.rho_dh = 0.f;
+  p->pressure_bar = 0.f;
+  p->density.pressure_bar_dh = 0.f;
 }
 
 /**
@@ -290,8 +307,13 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
  *
  * Multiplies the density and number of neighbours by the appropiate constants
  * and add the self-contribution term.
+ * Additional quantities such as velocity gradients will also get the final
+ * terms added to them here.
+ *
+ * Also adds/multiplies the cosmological terms if need be.
  *
  * @param p The particle to act upon
+ * @param cosmo The cosmological model.
  */
 __attribute__((always_inline)) INLINE static void hydro_end_density(
     struct part *restrict p, const struct cosmology *cosmo) {
@@ -301,40 +323,34 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   const float h_inv = 1.0f / h;                       /* 1/h */
   const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
   const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
-  
+
   /* Final operation on the density (add self-contribution). */
-  p->density.wcount += kernel_root;
-  p->density.wcount_dh -= hydro_dimension * kernel_root * h_inv;
   p->rho += p->mass * kernel_root;
+  p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
   p->pressure_bar += p->mass * p->u * kernel_root;
-  p->pressure_bar_dh -=
-    p->mass * p->u * h_inv * (hydro_dimension * kernel_root);
+  p->density.pressure_bar_dh -= hydro_dimension * p->mass * p->u * kernel_root;
+  p->density.wcount += kernel_root;
+  p->density.wcount_dh -= hydro_dimension * kernel_root;
 
   /* Finish the calculation by inserting the missing h-factors */
+  p->rho *= h_inv_dim;
+  p->density.rho_dh *= h_inv_dim_plus_one;
+  p->pressure_bar *= (h_inv_dim * hydro_gamma_minus_one);
+  p->density.pressure_bar_dh *= (h_inv_dim_plus_one * hydro_gamma_minus_one);
   p->density.wcount *= h_inv_dim;
-  p->density.wcount_dh *= h_inv_dim;
-  p->pressure_bar *= h_inv_dim;
-  p->pressure_bar_dh *= h_inv_dim;
-
-  /* Insert gamma factors */
-  p->pressure_bar *= hydro_gamma_minus_one;
-  p->pressure_bar_dh *= hydro_gamma_minus_one;
-
-  /* Finish calculation of the velocity curl components */
-  const float factor = h_inv_dim_plus_one / p->rho;
-  p->density.rot_v[0] *= factor;
-  p->density.rot_v[1] *= factor;
-  p->density.rot_v[2] *= factor;
-
-  /* Finish calculation of the velocity divergence */
-  p->density.div_v *= factor;
+  p->density.wcount_dh *= h_inv_dim_plus_one;
 }
 
 /**
  * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
  *
+ * In the desperate case where a particle has no neighbours (likely because
+ * of the h_max ceiling), set the particle fields to something sensible to avoid
+ * NaNs in the next calculations.
+ *
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
+ * @param cosmo The cosmological model.
  */
 __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
     struct part *restrict p, struct xpart *restrict xp,
@@ -346,68 +362,54 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
   const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
 
   /* Re-set problematic values */
+  p->rho = p->mass * kernel_root * h_inv_dim;
+  p->pressure_bar = p->mass * p->u * hydro_gamma_minus_one * kernel_root * h_inv_dim;
   p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+  p->density.rho_dh = 0.f;
   p->density.wcount_dh = 0.f;
-  p->pressure_bar =
-    kernel_root * p->mass * p->u * hydro_gamma_minus_one * h_inv_dim;
-  p->rho = kernel_root * p->mass * h_inv_dim;
-  p->pressure_bar_dh = 0;
-  p->density.div_v = 0.f;
-  p->density.rot_v[0] = 0.f;
-  p->density.rot_v[1] = 0.f;
-  p->density.rot_v[2] = 0.f;
+  p->density.pressure_bar_dh = 0.f;
 }
 
 /**
  * @brief Prepare a particle for the force calculation.
  *
- * Computes viscosity term, conduction term and smoothing length gradient terms.
+ * This function is called in the ghost task to convert some quantities coming
+ * from the density loop over neighbours into quantities ready to be used in the
+ * force loop over neighbours. Quantities are typically read from the density
+ * sub-structure and written to the force sub-structure.
+ * Examples of calculations done here include the calculation of viscosity term
+ * constants, thermal conduction terms, hydro conversions, etc.
  *
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
+ * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part *restrict p, struct xpart *restrict xp,
     const struct cosmology *cosmo) {
-  
-  const float fac_mu = 1.f; /* Will change with cosmological integration */
-
-  /* Compute the norm of the curl */
-  const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
-                             p->density.rot_v[1] * p->density.rot_v[1] +
-                             p->density.rot_v[2] * p->density.rot_v[2]);
-
-  /* Compute the norm of div v */
-  const float abs_div_v = fabsf(p->density.div_v);
 
   /* Compute the pressure */
-  const float pressure =
-    gas_pressure_from_internal_energy(p->rho, p->u);
+  const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
 
   /* Compute the sound speed */
   const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
 
-  /* Compute this `useful' property */
+  /* Compute the "grad h" term */
   const float rho_inv = 1.f / p->rho;
-  const float P_over_rho2 = pressure * rho_inv * rho_inv;
-
-  /* Compute the Balsara switch */
-  const float balsara =
-    abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);
-
+  const float grad_h_term =
+      1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
 
-  /* Update variables */
+  /* Update variables. */
+  p->force.f = grad_h_term;
+  p->force.pressure = pressure;
   p->force.soundspeed = soundspeed;
-  p->force.balsara = balsara;
-  p->force.P_over_rho2 = P_over_rho2;
-  p->force.f = 1.f; /* This cannot be computed individually in Pressure-Energy */
 }
 
 /**
  * @brief Reset acceleration fields of a particle
  *
  * Resets all hydro acceleration and time derivative fields in preparation
- * for the sums taking place in the variaous force tasks
+ * for the sums taking  place in the various force tasks.
  *
  * @param p The particle to act upon
  */
@@ -420,10 +422,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   p->a_hydro[2] = 0.0f;
 
   /* Reset the time derivatives. */
-  p->force.h_dt = 0.0f;
   p->u_dt = 0.0f;
-
-  /* Reset maximal signal velocity */
+  p->force.h_dt = 0.0f;
   p->force.v_sig = 0.0f;
 }
 
@@ -442,15 +442,23 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
   p->v[1] = xp->v_full[1];
   p->v[2] = xp->v_full[2];
 
+  /* Re-set the entropy */
   p->u = xp->u_full;
 }
 
 /**
  * @brief Predict additional particle fields forward in time when drifting
  *
- * @param p The particle
- * @param xp The extended data of the particle
- * @param dt The drift time-step.
+ * Additional hydrodynamic quantites are drifted forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * Note the different time-step sizes used for the different quantities as they
+ * include cosmological factors.
+ *
+ * @param p The particle.
+ * @param xp The extended data of the particle.
+ * @param dt_drift The drift time-step for positions.
+ * @param dt_therm The drift time-step for thermal quantities.
  */
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     struct part *restrict p, const struct xpart *restrict xp, float dt_drift,
@@ -466,99 +474,102 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     p->h *= expf(w1);
 
   /* Predict density */
-
   const float w2 = -hydro_dimension * w1;
-  if (fabsf(w2) < 0.2f) {
+  if (fabsf(w2) < 0.2f)
     p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
-    p->pressure_bar *= approx_expf(w2);
-  } else {
+  else
     p->rho *= expf(w2);
-    p->pressure_bar *= expf(w2);
-  }
 
-  /* Predict the energy */
+  /* Predict the internal energy */
   p->u += p->u_dt * dt_therm;
 
-  /* Compute the pressure */
+  /* Compute the new pressure */
   const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
 
   /* Compute the new sound speed */
   const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
 
+  p->force.pressure = pressure;
   p->force.soundspeed = soundspeed;
 }
 
 /**
  * @brief Finishes the force calculation.
  *
- * Multiplies the forces and accelerationsby the appropiate constants
+ * Multiplies the force and accelerations by the appropiate constants
+ * and add the self-contribution term. In most cases, there is little
+ * to do here.
+ *
+ * Cosmological terms are also added/multiplied here.
  *
  * @param p The particle to act upon
+ * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void hydro_end_force(
     struct part *restrict p, const struct cosmology *cosmo) {
-  
+
   p->force.h_dt *= p->h * hydro_dimension_inv;
 }
 
 /**
  * @brief Kick the additional variables
  *
- * @param p The particle to act upon
- * @param xp The particle extended data to act upon
- * @param dt The time-step for this kick
- * @param half_dt The half time-step for this kick
+ * Additional hydrodynamic quantites are kicked forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * @param p The particle to act upon.
+ * @param xp The particle extended data to act upon.
+ * @param dt_therm The time-step for this kick (for thermodynamic quantities).
  */
 __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     struct part *restrict p, struct xpart *restrict xp, float dt_therm) {
 
-  /* Do not decrease the energy (temperature) by more than a factor of 2*/
-  if (dt_therm > 0. && p->u_dt * dt_therm < -0.5f * p->u) {
-    p->u_dt = -0.5f * p->u / dt_therm;
+  /* Do not decrease the energy by more than a factor of 2*/
+  if (dt_therm > 0. && p->u_dt * dt_therm < -0.5f * xp->u_full) {
+    p->u_dt = -0.5f * xp->u_full / dt_therm;
   }
   xp->u_full += p->u_dt * dt_therm;
 
   /* Compute the pressure */
-  const float pressure =
-      gas_pressure_from_internal_energy(p->rho, xp->u_full);
+  const float pressure = gas_pressure_from_internal_energy(p->rho, xp->u_full);
 
-  /* Compute the new sound speed */
-  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+  /* Compute the sound speed */
+  const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
 
-  /* Update the sound speed */
+  p->force.pressure = pressure;
   p->force.soundspeed = soundspeed;
 }
 
 /**
- *  @brief Converts hydro quantity of a particle at the start of a run
+ * @brief Converts hydro quantity of a particle at the start of a run
  *
- * Requires the density to be known
+ * This function is called once at the end of the engine_init_particle()
+ * routine (at the start of a calculation) after the densities of
+ * particles have been computed.
+ * This can be used to convert internal energy into entropy for instance.
  *
  * @param p The particle to act upon
+ * @param xp The extended particle to act upon
  */
 __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
     struct part *restrict p, struct xpart *restrict xp) {
 
-  p->entropy = hydro_get_comoving_entropy(p);
-
   /* Compute the pressure */
   const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
-  /* Compute the sound speed */
-  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
 
-  const float rho_inv = 1.f / p->rho;
-  const float P_over_rho2 = pressure * rho_inv * rho_inv;
+  /* Compute the sound speed */
+  const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
 
-  /* Update sound speed */
+  p->force.pressure = pressure;
   p->force.soundspeed = soundspeed;
-  p->force.P_over_rho2 = P_over_rho2;
 }
 
 /**
  * @brief Initialises the particles for the first time
  *
  * This function is called only once just after the ICs have been
- * read in to do some conversions.
+ * read in to do some conversions or assignments between the particle
+ * and extended particle fields.
  *
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
@@ -567,25 +578,18 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
     struct part *restrict p, struct xpart *restrict xp) {
 
   p->time_bin = 0;
-  p->rho = 0.f;
-  p->pressure_bar = 0.f;
   xp->v_full[0] = p->v[0];
   xp->v_full[1] = p->v[1];
   xp->v_full[2] = p->v[2];
-
+  xp->a_grav[0] = 0.f;
+  xp->a_grav[1] = 0.f;
+  xp->a_grav[2] = 0.f;
   xp->u_full = p->u;
 
   hydro_reset_acceleration(p);
   hydro_init_part(p, NULL);
 }
 
-/* Begin actual helper functions for the specific hydro implementation.
- *
- * Here we define:
- *
- * hydro_h_term(*pi, *pj)
- */
-
 /**
  *  @brief Calculates the f_ij factors for a given particle i and j (i.e. the
  *         'h-terms'.
@@ -602,7 +606,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
  *  @param hi The smoothing length of particle i
  */
 __attribute__((always_inline)) INLINE static float hydro_h_term(
-    struct part *restrict pi, struct part *restrict pj, float hi) {
+    const struct part * pi, const struct part * pj, float hi) {
   
   /* First constrct h_i / (n_D n_i) */
   
@@ -612,7 +616,7 @@ __attribute__((always_inline)) INLINE static float hydro_h_term(
   const float bottom_term = 1 + hi_over_density * pi->density.wcount_dh;
 
   /* Construct the 'top term' */
-  const float top_term_top = pi->pressure_bar_dh * hi_over_density;
+  const float top_term_top = pi->density.pressure_bar_dh * hi_over_density;
   const float top_term_bottom =
     hydro_gamma_minus_one * pj->mass * pj->u;
 
@@ -622,4 +626,5 @@ __attribute__((always_inline)) INLINE static float hydro_h_term(
   return f_ij;
 }
 
-#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_H */
+
+#endif /* SWIFT_MINIMAL_HYDRO_H */
diff --git a/src/hydro/PressureEnergy/hydro_debug.h b/src/hydro/PressureEnergy/hydro_debug.h
index 098954f3ad9619a143f391a904b2a5de461e3e7b..41d295a5b76b7ebc4b297271fe908b38428bdd16 100644
--- a/src/hydro/PressureEnergy/hydro_debug.h
+++ b/src/hydro/PressureEnergy/hydro_debug.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published
@@ -18,29 +18,35 @@
  ******************************************************************************/
 #ifndef SWIFT_PRESSURE_ENERGY_HYDRO_DEBUG_H
 #define SWIFT_PRESSURE_ENERGY_HYDRO_DEBUG_H
-
 /**
- * @file PressureEnergy/hydro_debug.h
- * @brief Pressure-Energy implementation of SPH (Debugging routines)
+ * @file Minimal/hydro_debug.h
+ * @brief Minimal conservative implementation of SPH (Debugging routines)
  *
- * The thermal variable is the energy (u) and the pressure is smoothed over
- * contact discontinuities to prevent spurious surface tension.
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
  *
- * Follows equations (16), (17) and (18) of Hopkins, P., MNRAS, 2013,
- * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
+ * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
  */
 
 __attribute__((always_inline)) INLINE static void hydro_debug_particle(
     const struct part* p, const struct xpart* xp) {
   printf(
       "x=[%.3e,%.3e,%.3e], "
-      "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
-      "h=%.3e, wcount=%.3f, wcount_dh=%.3e, time_bin=%d\n"
-      "p_bar=%.3e, p_bar_dh=%3.e, u=%3.e\n",
+      "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
+      "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
+      "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, \n"
+      "p_dh=%.3e, p_bar=%.3e \n"
+      "time_bin=%d\n",
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
       xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
-      p->h, p->density.wcount, p->density.wcount_dh, p->time_bin,
-      p->pressure_bar, p->pressure_bar_dh, p->u);
+      p->u, p->u_dt, p->force.v_sig, hydro_get_comoving_pressure(p), p->h,
+      p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
+      p->density.pressure_bar_dh, p->pressure_bar,
+      p->time_bin);
 }
 
-#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_DEBUG_H */
+#endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */
diff --git a/src/hydro/PressureEnergy/hydro_iact.h b/src/hydro/PressureEnergy/hydro_iact.h
index 7ee84c8128f14bce45910cebc87790a404e96c42..70ba6e784823b7892d87e3682adf36b7b0381bdb 100644
--- a/src/hydro/PressureEnergy/hydro_iact.h
+++ b/src/hydro/PressureEnergy/hydro_iact.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * This file is part* of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published
@@ -16,201 +16,323 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H
-#define SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H
+#ifndef SWIFT_MINIMAL_HYDRO_IACT_H
+#define SWIFT_MINIMAL_HYDRO_IACT_H
 
 /**
- * @file PressureEnergy/hydro_iact.h
- * @brief Pressure-Energy implementation of SPH (Neighbour loop equations)
+ * @file Minimal/hydro_iact.h
+ * @brief Minimal conservative implementation of SPH (Neighbour loop equations)
  *
- * The thermal variable is the energy (u) and the pressure is smoothed over
- * contact discontinuities to prevent spurious surface tension.
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
  *
- * Follows equations (16), (17) and (18) of Hopkins, P., MNRAS, 2013,
- * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
+ * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational
+ * Physics, 2012, Volume 231, Issue 3, pp. 759-794.
  */
 
+#include "adiabatic_index.h"
+#include "minmax.h"
+
 /**
- * @brief Density loop (non-symmetric version)
+ * @brief Density interaction between two part*icles.
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
  */
-__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
-    float a, float H) {
+__attribute__((always_inline)) INLINE static void runner_iact_density(
+    float r2, const float *dx, float hi, float hj, struct part* pi,
+    struct part* pj, float a, float H) {
 
-  float wi, wi_dx;
-  float dv[3], curlvr[3]; 
+  float wi, wj, wi_dx, wj_dx;
 
-  /* Get r and r inverse. */
-  const float r_inv = 1.0f / sqrtf(r2);
-  const float r = r2 * r_inv;
+  const float r = sqrtf(r2);
 
-  /* Get the masses and energies */
+  /* Get the masses. */
+  const float mi = pi->mass;
   const float mj = pj->mass;
-  const float int_energy_j = pj->u;
 
-  /* Compute the kernel function */
-  const float hi_inv = 1.0f / hi;
-  const float hj_inv = 1.0f / hj;
+  /* Compute density of pi. */
+  const float hi_inv = 1.f / hi;
   const float ui = r * hi_inv;
   kernel_deval(ui, &wi, &wi_dx);
 
-  /* Compute contribution to the number of neighbours */
-  /* Note that wcount is \bar{n} */
+  pi->rho += mj * wi;
+  pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+  pi->pressure_bar += mj * wi * pj->u;
+  pi->density.pressure_bar_dh -= mj * pj->u * (hydro_dimension * wi + ui * wi_dx);
   pi->density.wcount += wi;
-  pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx) * hj_inv;
-
-  /* Compute the 'real' SPH density for other uses */
-  const float wimj = wi * mj;
-  pi->rho += wimj;
-
-  /* Compute the contribution to the local pressure */
-  /* The (gamma - 1) factor is added in hydro_end_density. */
-  pi->pressure_bar += int_energy_j * wimj;
-  pi->pressure_bar_dh -=
-    mj * int_energy_j * hj_inv * (hydro_dimension * wi + ui * wi_dx);
-
-  /* The following is lifted directly from the PressureEntropy code */
-  const float fac = mj * wi_dx * r_inv;
-
-  /* Compute dv dot r */
-  dv[0] = pi->v[0] - pj->v[0];
-  dv[1] = pi->v[1] - pj->v[1];
-  dv[2] = pi->v[2] - pj->v[2];
-  const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
-  pi->density.div_v -= fac * dvdr;
-
-  /* Compute dv cross r */
-  curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
-  curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
-  curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
-
-  pi->density.rot_v[0] += fac * curlvr[0];
-  pi->density.rot_v[1] += fac * curlvr[1];
-  pi->density.rot_v[2] += fac * curlvr[2];
+  pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+  /* Compute density of pj. */
+  const float hj_inv = 1.f / hj;
+  const float uj = r * hj_inv;
+  kernel_deval(uj, &wj, &wj_dx);
+
+  pj->rho += mi * wj;
+  pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
+  pj->pressure_bar += mi * wj * pi->u;
+  pj->density.pressure_bar_dh -= mi * pi->u * (hydro_dimension * wj + uj * wj_dx);
+  pj->density.wcount += wj;
+  pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
 }
 
 /**
- * @brief Density loop
+ * @brief Density interaction between two part*icles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
  */
-__attribute__((always_inline)) INLINE static void runner_iact_density(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
-    float a, float H) {
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
+    float r2, const float *dx, float hi, float hj, struct part* pi,
+    const struct part* pj, float a, float H) {
+
+  float wi, wi_dx;
+
+  /* Get the masses. */
+  const float mj = pj->mass;
 
-  /* For now, just do the simple case */
-  float negative_dx[3];
+  /* Get r and r inverse. */
+  const float r = sqrtf(r2);
 
-  negative_dx[0] = -1.f * dx[0];
-  negative_dx[1] = -1.f * dx[1];
-  negative_dx[2] = -1.f * dx[2];
+  const float h_inv = 1.f / hi;
+  const float ui = r * h_inv;
+  kernel_deval(ui, &wi, &wi_dx);
 
-  runner_iact_nonsym_density(r2, dx, hi, hj, pi, pj, a, H);
-  runner_iact_nonsym_density(r2, negative_dx, hj, hi, pj, pi, a, H);
+  pi->rho += mj * wi;
+  pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+  pi->pressure_bar += mj * wi * pj->u;
+  pi->density.pressure_bar_dh -= mj * pj->u * (hydro_dimension * wi + ui * wi_dx);
+  pi->density.wcount += wi;
+  pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
 }
 
 /**
- * @brief Force loop (non-symmetric version)
+ * @brief Force interaction between two part*icles.
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
  */
-__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
-    float a, float H) {
-  
-  float wi, wi_dx, wj, wj_dx;
-  const float fac_mu = 1.f; /* Will change with cosmological integration */
-  float dv[3];
-
-  /* Masses */
-  const float mj = pj->mass;
+__attribute__((always_inline)) INLINE static void runner_iact_force(
+    float r2, const float *dx, float hi, float hj, struct part* pi,
+    struct part* pj, float a, float H) {
 
-  /* Get r and r inverse */
-  const float r = sqrtf(r2);
-  const float ri = 1.0f / r;
+  /* Cosmological factors entering the EoMs */
+  const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+  const float a2_Hubble = a * a * H;
 
-  /* Get velocity difference */
-  dv[0] = pi->v[0] - pj->v[0];
-  dv[1] = pi->v[1] - pj->v[1];
-  dv[2] = pi->v[2] - pj->v[2];
+  const float r = sqrtf(r2);
+  const float r_inv = 1.0f / r;
 
-  const float dvdr = dv[0] * dx[0] +
-                     dv[1] * dx[1] +
-                     dv[2] * dx[2];
+  /* Recover some data */
+  const float mi = pi->mass;
+  const float mj = pj->mass;
+  const float rhoi = pi->rho;
+  const float rhoj = pj->rho;
+  /* Compute f terms */
+  const float f_ij = hydro_h_term(pi, pj, hi);
+  const float f_ji = hydro_h_term(pj, pi, hj);
 
-  /* Compute kernel function */
+  /* Get the kernel for hi. */
   const float hi_inv = 1.0f / hi;
+  const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+  const float xi = r * hi_inv;
+  float wi, wi_dx;
+  kernel_deval(xi, &wi, &wi_dx);
+  const float wi_dr = hid_inv * wi_dx;
+
+  /* Get the kernel for hj. */
   const float hj_inv = 1.0f / hj;
-  const float hid_inv = pow_dimension_plus_one(hi_inv);
-  const float hjd_inv = pow_dimension_plus_one(hj_inv);
+  const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
+  const float xj = r * hj_inv;
+  float wj, wj_dx;
+  kernel_deval(xj, &wj, &wj_dx);
+  const float wj_dr = hjd_inv * wj_dx;
 
-  const float ui = r * hi_inv;
-  const float uj = r * hj_inv;
-  kernel_deval(ui, &wi, &wi_dx);
-  kernel_deval(uj, &wj, &wj_dx);
+  /* Compute dv dot r. */
+  const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+                     (pi->v[1] - pj->v[1]) * dx[1] +
+                     (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
 
-  const float wi_dr = hid_inv * wi_dx;
-  const float wj_dr = hjd_inv * wj_dx;
+  /* Are the part*icles moving towards each others ? */
+  const float omega_ij = min(dvdr, 0.f);
+  const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
 
-  /* Compute gradient terms */
+  /* Compute sound speeds and signal velocity */
+  const float ci = pi->force.soundspeed;
+  const float cj = pj->force.soundspeed;
+  const float v_sig = ci + cj - 3.f * mu_ij;
+
+  /* Construct the full viscosity term */
+  const float rho_ij = 0.5f * (rhoi + rhoj);
+  const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+
+  /* Convolve with the kernel */
+  const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+
+  /* SPH acceleration term */
+  const float sph_acc_term = pj->u * pi->u *
+    hydro_gamma_minus_one * hydro_gamma_minus_one *
+    ((f_ij/pi->pressure_bar) * wi_dr + (f_ji/pj->pressure_bar) * wj_dr) * r_inv; 
+
+  /* Assemble the acceleration */
+  const float acc = sph_acc_term + visc_acc_term;
+
+  /* Use the force Luke ! */
+  pi->a_hydro[0] -= mj * acc * dx[0];
+  pi->a_hydro[1] -= mj * acc * dx[1];
+  pi->a_hydro[2] -= mj * acc * dx[2];
+
+  pj->a_hydro[0] += mi * acc * dx[0];
+  pj->a_hydro[1] += mi * acc * dx[1];
+  pj->a_hydro[2] += mi * acc * dx[2];
+
+  /* Get the time derivative for u. */
+  const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
+    pj->u * pi->u * (f_ij / pi->pressure_bar) * wi_dr * dvdr * r_inv;
+  const float sph_du_term_j = hydro_gamma_minus_one * hydro_gamma_minus_one *
+    pi->u * pj->u * (f_ji / pj->pressure_bar) * wj_dr * dvdr * r_inv;
+
+  /* Viscosity term */
+  const float visc_du_term = 0.5f * visc_acc_term * dvdr;
+
+  /* Assemble the energy equation term */
+  const float du_dt_i = sph_du_term_i + visc_du_term;
+  const float du_dt_j = sph_du_term_j + visc_du_term;
+
+  /* Internal energy time derivatibe */
+  pi->u_dt += du_dt_i * mj;
+  pj->u_dt += du_dt_j * mi;
+
+  /* Get the time derivative for h. */
+  pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+  pj->force.h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
+
+  /* Update the signal velocity. */
+  pi->force.v_sig = max(pi->force.v_sig, v_sig);
+  pj->force.v_sig = max(pj->force.v_sig, v_sig);
+}
+
+/**
+ * @brief Force interaction between two part*icles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
+    float r2, const float *dx, float hi, float hj, struct part* pi,
+    const struct part* pj, float a, float H) {
 
+  /* Cosmological factors entering the EoMs */
+  const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+  const float a2_Hubble = a * a * H;
+
+  const float r = sqrtf(r2);
+  const float r_inv = 1.0f / r;
+
+  /* Recover some data */
+  // const float mi = pi->mass;
+  const float mj = pj->mass;
+  const float rhoi = pi->rho;
+  const float rhoj = pj->rho;
+  /* Compute gradient terms */
   const float f_ij = hydro_h_term(pi, pj, hi);
   const float f_ji = hydro_h_term(pj, pi, hj);
 
-  /* Artificial viscosity */
 
-  const float omega_ij = fminf(dvdr, 0.f);
-  const float mu_ij = fac_mu * ri * omega_ij;
+  /* Get the kernel for hi. */
+  const float hi_inv = 1.0f / hi;
+  const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+  const float xi = r * hi_inv;
+  float wi, wi_dx;
+  kernel_deval(xi, &wi, &wi_dx);
+  const float wi_dr = hid_inv * wi_dx;
 
-  const float v_sig = pi->force.soundspeed + pj->force.soundspeed - 3.f * mu_ij;
+  /* Get the kernel for hj. */
+  const float hj_inv = 1.0f / hj;
+  const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
+  const float xj = r * hj_inv;
+  float wj, wj_dx;
+  kernel_deval(xj, &wj, &wj_dx);
+  const float wj_dr = hjd_inv * wj_dx;
 
-  const float rho_ij = 0.5f * (pi->rho + pj->rho);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (pi->force.balsara + pj->force.balsara) / rho_ij;
-  const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
+  /* Compute dv dot r. */
+  const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+                     (pi->v[1] - pj->v[1]) * dx[1] +
+                     (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
 
-  /* Calculate the equation of motion H13:15 */
+  /* Are the part*icles moving towards each others ? */
+  const float omega_ij = min(dvdr, 0.f);
+  const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
 
-  const float sph_term = mj * pj->u * pi->u *
-    hydro_gamma_minus_one * hydro_gamma_minus_one *
-    ((f_ij/pi->pressure_bar) * wi_dr + (f_ji/pj->pressure_bar) * wj_dr);
+  /* Compute sound speeds and signal velocity */
+  const float ci = pi->force.soundspeed;
+  const float cj = pj->force.soundspeed;
+  const float v_sig = ci + cj - 3.f * mu_ij;
 
-  const float acc = (visc_term + sph_term) * ri;
+  /* Construct the full viscosity term */
+  const float rho_ij = 0.5f * (rhoi + rhoj);
+  const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
 
-  /* Use the force Luke ! */
-  pi->a_hydro[0] -= acc * dx[0];
-  pi->a_hydro[1] -= acc * dx[1];
-  pi->a_hydro[2] -= acc * dx[2];
+  /* Convolve with the kernel */
+  const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
 
-  /* Get the time derivative for h */
-  pi->force.h_dt -= mj * dvdr * ri / pj->rho * wi_dr;
+  /* SPH acceleration term */
+  const float sph_acc_term = pj->u * pi->u *
+    hydro_gamma_minus_one * hydro_gamma_minus_one *
+    ((f_ij/pi->pressure_bar) * wi_dr + (f_ji/pj->pressure_bar) * wj_dr) * r_inv; 
 
-  /* Update the signal velocity */
-  pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
+  /* Assemble the acceleration */
+  const float acc = sph_acc_term + visc_acc_term;
 
-  /* Calcualte change in energy from viscosity */
-  const float u_dt_visc = mj * visc_term * dvdr * ri;
+  /* Use the force Luke ! */
+  pi->a_hydro[0] -= mj * acc * dx[0];
+  pi->a_hydro[1] -= mj * acc * dx[1];
+  pi->a_hydro[2] -= mj * acc * dx[2];
 
-  /* Calculate the change in internal energy from hydro H13:16 */
-  const float u_dt_hydro =
-    hydro_gamma_minus_one * hydro_gamma_minus_one *
-    mj * pj->u * pi->u * (f_ij / pi->pressure_bar) * wi_dr * dvdr * ri;
+  /* Get the time derivative for u. */
+  const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
+    pj->u * pi->u * (f_ij / pi->pressure_bar) * wi_dr * dvdr * r_inv;
 
-  pi->u_dt += (u_dt_hydro + u_dt_visc);
-}
+  /* Viscosity term */
+  const float visc_du_term = 0.5f * visc_acc_term * dvdr;
 
-/**
- * @brief Force loop
- */
-__attribute__((always_inline)) INLINE static void runner_iact_force(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
-    float a, float H) {
+  /* Assemble the energy equation term */
+  const float du_dt_i = sph_du_term_i + visc_du_term;
 
-  /* For now, just do the simple case */
-  float negative_dx[3];
+  /* Internal energy time derivatibe */
+  pi->u_dt += du_dt_i * mj;
 
-  negative_dx[0] = -1.f * dx[0];
-  negative_dx[1] = -1.f * dx[1];
-  negative_dx[2] = -1.f * dx[2];
+  /* Get the time derivative for h. */
+  pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
 
-  runner_iact_nonsym_force(r2, dx, hi, hj, pi, pj, a, H);
-  runner_iact_nonsym_force(r2, negative_dx, hj, hi, pj, pi, a, H);
+  /* Update the signal velocity. */
+  pi->force.v_sig = max(pi->force.v_sig, v_sig);
 }
 
-#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H */
+#endif /* SWIFT_MINIMAL_HYDRO_IACT_H */
diff --git a/src/hydro/PressureEnergy/hydro_io.h b/src/hydro/PressureEnergy/hydro_io.h
index 530f387f17063f9a033b6bbff25cca10afd1ca88..0bc0534172c5bb0764cd0a2ddb2024d2cb9eb6a0 100644
--- a/src/hydro/PressureEnergy/hydro_io.h
+++ b/src/hydro/PressureEnergy/hydro_io.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published
@@ -16,18 +16,20 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef SWIFT_PRESSURE_ENERGY_HYDRO_IO_H
-#define SWIFT_PRESSURE_ENERGY_HYDRO_IO_H
-
+#ifndef SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H
+#define SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H
 /**
- * @file PressureEnergy/hydro_io.h
- * @brief Pressure-Energy implementation of SPH (i/o routines)
+ * @file Minimal/hydro_io.h
+ * @brief Minimal conservative implementation of SPH (i/o routines)
  *
- * The thermal variable is the energy (u) and the pressure is smoothed over
- * contact discontinuities to prevent spurious surface tension.
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
  *
- * Follows equations (17) and (18) of Hopkins, P., MNRAS, 2013,
- * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
+ * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
  */
 
 #include "adiabatic_index.h"
@@ -52,21 +54,30 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
                                 UNIT_CONV_LENGTH, parts, x);
   list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
                                 UNIT_CONV_SPEED, parts, v);
-  list[2] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
+  list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+                                parts, mass);
+  list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
                                 UNIT_CONV_LENGTH, parts, h);
-  list[3] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
-                                 parts, mass);
-  list[4] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+  list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
+                                UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
+  list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
                                 UNIT_CONV_NO_UNITS, parts, id);
-  list[5] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
-                                UNIT_CONV_ENERGY_PER_UNIT_MASS, parts,
-                                u);
   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, rho);
 }
 
+void convert_S(const struct engine* e, const struct part* p, float* ret) {
+
+  ret[0] = hydro_get_comoving_entropy(p);
+}
+
+void convert_P(const struct engine* e, const struct part* p, float* ret) {
+
+  ret[0] = hydro_get_comoving_pressure(p);
+}
+
 void convert_part_pos(const struct engine* e, const struct part* p,
                       double* ret) {
 
@@ -94,26 +105,24 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
   *num_fields = 9;
 
   /* List what we want to write */
-  list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
-                                              UNIT_CONV_LENGTH, parts,
-                                              convert_part_pos);
-  list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED,
-                                 parts, v);
-  list[2] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+  list[0] = io_make_output_field_convert_part(
+      "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
+  list[1] =
+      io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
+  list[2] =
+      io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
+  list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
                                  parts, h);
-  list[3] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
-                                 UNIT_CONV_NO_UNITS, parts, id);
   list[4] = io_make_output_field("InternalEnergy", FLOAT, 1,
-                                 UNIT_CONV_ENERGY_PER_UNIT_MASS,
-                                 parts, u);
-  list[5] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts,
-                                 mass);
-  list[6] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
-                                 rho);
-  list[7] = io_make_output_field("Accelerations", FLOAT, 3,
-                                 UNIT_CONV_ACCELERATION, parts, a_hydro);
-  list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
-                                 parts, pressure_bar);
+                                 UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
+  list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
+                                 UNIT_CONV_NO_UNITS, parts, id);
+  list[6] =
+      io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
+  list[7] = io_make_output_field_convert_part(
+      "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, convert_S);
+  list[8] = io_make_output_field_convert_part(
+      "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
 }
 
 /**
@@ -121,12 +130,12 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
  * @param h_grpsph The HDF5 group in which to write
  */
 void hydro_write_flavour(hid_t h_grpsph) {
+
+  /* Viscosity and thermal conduction */
+  /* Nothing in this minimal model... */
   io_write_attribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
-  io_write_attribute_s(
-      h_grpsph, "Viscosity Model",
-      "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
-  io_write_attribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
-  io_write_attribute_f(h_grpsph, "Viscosity beta", 3.f);
+  io_write_attribute_s(h_grpsph, "Viscosity Model",
+                       "Minimal treatment as in Monaghan (1992)");
 
   /* Time integration properties */
   io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
@@ -140,4 +149,4 @@ void hydro_write_flavour(hid_t h_grpsph) {
  */
 int writeEntropyFlag() { return 0; }
 
-#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_IO_H */
+#endif /* SWIFT_MINIMAL_HYDRO_IO_H */
diff --git a/src/hydro/PressureEnergy/hydro_part.h b/src/hydro/PressureEnergy/hydro_part.h
index 65a927873ad2a3b526c8069be56e947e51470c9a..cd09b7d876f58d5bc96bc39cf910c5a87453f403 100644
--- a/src/hydro/PressureEnergy/hydro_part.h
+++ b/src/hydro/PressureEnergy/hydro_part.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published
@@ -18,37 +18,44 @@
  ******************************************************************************/
 #ifndef SWIFT_PRESSURE_ENERGY_HYDRO_PART_H
 #define SWIFT_PRESSURE_ENERGY_HYDRO_PART_H
-
 /**
- * @file PressureEnergy/hydro_part.h
- * @brief Pressure-Energy implementation of SPH (Particle definition)
+ * @file Minimal/hydro_part.h
+ * @brief Minimal conservative implementation of SPH (Particle definition)
  *
- * The thermal variable is the energy (u) and the pressure is smoothed over
- * contact discontinuities to prevent spurious surface tension.
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
  *
- * Follows equations (16), (17) and (18) of Hopkins, P., MNRAS, 2013,
- * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
+ * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational
+ * Physics, 2012, Volume 231, Issue 3, pp. 759-794.
  */
 
 #include "chemistry_struct.h"
 #include "cooling_struct.h"
 
-/* Extra particle data not needed during the SPH loops over neighbours. */
+/**
+ * @brief Particle fields not needed during the SPH loops over neighbours.
+ *
+ * This structure contains the particle fields that are not used in the
+ * density or force loops. Quantities should be used in the kick, drift and
+ * potentially ghost tasks only.
+ */
 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. */
+  /*! 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 */
+  /*! Gravitational acceleration at the last full step. */
   float a_grav[3];
 
-  /*! Full internal energy */
+  /*! Internal energy at the last full step. */
   float u_full;
 
   /*! Additional data used to record cooling information */
@@ -56,10 +63,16 @@ struct xpart {
 
 } SWIFT_STRUCT_ALIGN;
 
-/* Data of a single particle. */
+/**
+ * @brief Particle fields for the SPH particles
+ *
+ * The density and force substructures are used to contain variables only used
+ * within the density and force loops over neighbours. All more permanent
+ * variables should be declared in the main part of the part structure,
+ */
 struct part {
 
-  /*! Particle ID. */
+  /*! Particle unique ID. */
   long long id;
 
   /*! Pointer to corresponding gravity part. */
@@ -77,77 +90,78 @@ struct part {
   /*! Particle mass. */
   float mass;
 
-  /*! SPH Density */
-  float rho;
-
-  /*! Smoothed particle pressure. */
-  float pressure_bar;
-
-  /*! Smoothed particle pressure's spatial derivative */
-  float pressure_bar_dh;
+  /*! Particle smoothing length. */
+  float h;
 
-  /*! Particle internal energy */
+  /*! Particle internal energy. */
   float u;
 
-  /*! Differential of the internal energy with respect to time */
+  /*! Time derivative of the internal energy. */
   float u_dt;
 
-  /*! Entropy (for if people want it) */ 
-  float entropy;
-
-  /*! Particle cutoff radius. */
-  float h;
+  /*! Particle density. */
+  float rho;
+  
+  /*! Particle pressure (weighted) */
+  float pressure_bar;
 
+  /* Store density/force specific stuff. */
+  /*union {*/
 
+    /**
+     * @brief Structure for the variables only used in the density loop over
+     * neighbours.
+     *
+     * Quantities in this sub-structure should only be accessed in the density
+     * loop over neighbours and the ghost task.
+     */
     struct {
 
-      /*! Number of neighbours. */
+      /*! Neighbour number count. */
       float wcount;
 
-      /*! Number of neighbours spatial derivative. */
+      /*! Derivative of the neighbour number with respect to h. */
       float wcount_dh;
 
-      /*! Velocity curl */
-      float rot_v[3];
-
-      /*! Velocity divergence */
-      float div_v;
-
-      /*! d\rho/dh */
+      /*! Derivative of density with respect to h */
       float rho_dh;
 
-      /*! dP/dh */
-      float pressure_dh;
+	  /*! Derivative of the weighted pressure with respect to h */
+      float pressure_bar_dh;
 
     } density;
 
+    /**
+     * @brief Structure for the variables only used in the force loop over
+     * neighbours.
+     *
+     * Quantities in this sub-structure should only be accessed in the force
+     * loop over neighbours and the ghost, drift and kick tasks.
+     */
     struct {
 
-      /*! Signal velocity. */
-      float v_sig;
+      /*! "Grad h" term */
+      float f;
 
-      /*! Time derivative of the smoothing length */
-      float h_dt;
+      /*! Particle pressure. */
+      float pressure;
 
-      /*! Sound speed */
+      /*! Particle soundspeed. */
       float soundspeed;
-      
-      /*! Balsara switch */
-      float balsara;
 
-      /*! F_{ij} -- not actually possible to set this. */
-      float f;
+      /*! Particle signal velocity */
+      float v_sig;
 
-      /*! P/\rho^2 -- not actually required for Pressure Energy but this is
-      needed for cross-compatibility with the unit tests */
-      float P_over_rho2;
+      /*! Time derivative of smoothing length  */
+      float h_dt;
 
     } force;
+  /*};*/
 
   /* Chemistry information */
   struct chemistry_part_data chemistry_data;
 
-  /* Time-step length */
+  /*! Time-step length */
   timebin_t time_bin;
 
 #ifdef SWIFT_DEBUG_CHECKS
@@ -162,4 +176,4 @@ struct part {
 
 } SWIFT_STRUCT_ALIGN;
 
-#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_PART_H */
+#endif /* SWIFT_MINIMAL_HYDRO_PART_H */