diff --git a/configure.ac b/configure.ac
index f7b999354cde98dbda14b950cf2cea0c8a3f29fb..be7d82fa3bb71cec3511a6e5370061187afb80a5 100644
--- a/configure.ac
+++ b/configure.ac
@@ -425,6 +425,18 @@ elif test "$fixed_boundary_particles" != "no"; then
    AC_DEFINE_UNQUOTED([SWIFT_FIXED_BOUNDARY_PARTICLES], [$enableval] ,[Particles with smaller ID than this will be considered as boundaries.])
 fi
 
+# Check if fixed entropy is on for settling planetary initial conditions
+AC_ARG_ENABLE([planetary-fixed-entropy],
+   [AS_HELP_STRING([--enable-planetary-fixed-entropy],
+     [Force entropies to stay fixed for settling planetary initial conditions @<:@yes/no@:>@]
+   )],
+   [planetary_fixed_entropy="$enableval"],
+   [planetary_fixed_entropy="no"]
+)
+if test "$planetary_fixed_entropy" = "yes"; then
+   AC_DEFINE([PLANETARY_FIXED_ENTROPY],1,[Enable planetary fixed entropy])
+fi
+
 # Check whether we have any of the ARM v8.1 tick timers
 AX_ASM_ARM_PMCCNTR
 AX_ASM_ARM_CNTVCT
@@ -2468,6 +2480,7 @@ AC_MSG_RESULT([
    Custom icbrtf               : $enable_custom_icbrtf
    Boundary particles          : $boundary_particles
    Fixed boundary particles    : $fixed_boundary_particles
+   Planetary fixed entropy     : $planetary_fixed_entropy
 
    Particle Logger      : $with_logger
    Python enabled       : $have_python
diff --git a/doc/RTD/source/Planetary/index.rst b/doc/RTD/source/Planetary/index.rst
index 7b0a227e1bf387d383cfe42fa63d21c57ec87281..5cc45fa0e51fbfc38e42e7cf111a74ccfe9f835e 100644
--- a/doc/RTD/source/Planetary/index.rst
+++ b/doc/RTD/source/Planetary/index.rst
@@ -26,7 +26,7 @@ These allow for multiple materials to be used,
 chosen from the several available equations of state.
 
 .. toctree::
-   :maxdepth: 2
+   :maxdepth: 1
    :caption: More information:
    
    Hydro Scheme <hydro_scheme>
diff --git a/doc/RTD/source/Planetary/initial_conditions.rst b/doc/RTD/source/Planetary/initial_conditions.rst
index d4c20f1fa1074c57fa3864f6859980e5f1d832b0..fd348a9fd134fca3a35be4a7eeb16669ea56534a 100755
--- a/doc/RTD/source/Planetary/initial_conditions.rst
+++ b/doc/RTD/source/Planetary/initial_conditions.rst
@@ -20,4 +20,19 @@ Ruiz-Bonilla et al. (2020).
 They are available with documentation and examples at 
 https://github.com/srbonilla/WoMa and https://github.com/jkeger/seagen,
 or can be installed directly with ``pip``
-(https://pypi.org/project/woma/, https://pypi.org/project/seagen/).
\ No newline at end of file
+(https://pypi.org/project/woma/, https://pypi.org/project/seagen/).
+
+
+Settling initial conditions with fixed entropies
+------------------------------------------------
+
+If the particles' equations of state include specific entropies, 
+and the initial conditions file includes specific entropies for each particle
+(in ``PartType0/Entropies``), 
+then configuring SWIFT with ``--enable-planetary-fixed-entropy``
+will override the internal energy of each particle each step such that its 
+specific entropy remains constant. 
+
+This should be used with caution, but may be a convenient way to maintain an 
+entropy profile while initial conditions settle to equilibrium with their 
+slightly different SPH densities.
diff --git a/src/equation_of_state/planetary/sesame.h b/src/equation_of_state/planetary/sesame.h
index 2af28baebf333d77b5ccacafe1d8752624c8c93a..3480de806a1c665585bd6e603448175177e4be87 100644
--- a/src/equation_of_state/planetary/sesame.h
+++ b/src/equation_of_state/planetary/sesame.h
@@ -47,9 +47,9 @@ struct SESAME_params {
   float *table_log_u_rho_T;
   float *table_P_rho_T;
   float *table_c_rho_T;
-  float *table_s_rho_T;
+  float *table_log_s_rho_T;
   int date, num_rho, num_T;
-  float P_tiny, c_tiny;
+  float u_tiny, P_tiny, c_tiny, s_tiny;
   enum eos_planetary_material_id mat_id;
 };
 
@@ -138,7 +138,7 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat,
       (float *)malloc(mat->num_rho * mat->num_T * sizeof(float));
   mat->table_c_rho_T =
       (float *)malloc(mat->num_rho * mat->num_T * sizeof(float));
-  mat->table_s_rho_T =
+  mat->table_log_s_rho_T =
       (float *)malloc(mat->num_rho * mat->num_T * sizeof(float));
 
   // Densities (not log yet)
@@ -159,7 +159,8 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat,
     if (c != 1) error("Failed to read the SESAME EoS table %s", table_file);
   }
 
-  // Sp. int. energies (not log yet), pressures, sound speeds, and entropies
+  // Sp. int. energies (not log yet), pressures, sound speeds, and sp.
+  // entropies (not log yet)
   for (int i_T = -1; i_T < mat->num_T; i_T++) {
     for (int i_rho = -1; i_rho < mat->num_rho; i_rho++) {
       // Ignore the first elements of rho = 0, T = 0
@@ -171,7 +172,7 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat,
                    &mat->table_log_u_rho_T[i_rho * mat->num_T + i_T],
                    &mat->table_P_rho_T[i_rho * mat->num_T + i_T],
                    &mat->table_c_rho_T[i_rho * mat->num_T + i_T],
-                   &mat->table_s_rho_T[i_rho * mat->num_T + i_T]);
+                   &mat->table_log_s_rho_T[i_rho * mat->num_T + i_T]);
         if (c != 4) error("Failed to read the SESAME EoS table %s", table_file);
       }
     }
@@ -188,22 +189,11 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) {
     mat->table_log_rho[i_rho] = logf(mat->table_log_rho[i_rho]);
   }
 
-  // Convert sp. int. energies to log(sp. int. energy)
-  for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
-    for (int i_T = 0; i_T < mat->num_T; i_T++) {
-      // If not positive then set very small for the log
-      if (mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] <= 0) {
-        mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] = 1.f;
-      }
-
-      mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] =
-          logf(mat->table_log_u_rho_T[i_rho * mat->num_T + i_T]);
-    }
-  }
-
-  // Initialise tiny pressure and soundspeed
+  // Initialise tiny values
+  mat->u_tiny = FLT_MAX;
   mat->P_tiny = FLT_MAX;
   mat->c_tiny = FLT_MAX;
+  mat->s_tiny = FLT_MAX;
 
   // Enforce that the 1D arrays of u (at each rho) are monotonic
   // This is necessary because, for some high-density u slices at very low T,
@@ -223,7 +213,11 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) {
         break;
       }
 
-      // Smallest positive pressure and sound speed
+      // Smallest positive values
+      if ((mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] < mat->u_tiny) &&
+          (mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] > 0)) {
+        mat->u_tiny = mat->table_log_u_rho_T[i_rho * mat->num_T + i_T];
+      }
       if ((mat->table_P_rho_T[i_rho * mat->num_T + i_T] < mat->P_tiny) &&
           (mat->table_P_rho_T[i_rho * mat->num_T + i_T] > 0)) {
         mat->P_tiny = mat->table_P_rho_T[i_rho * mat->num_T + i_T];
@@ -232,12 +226,38 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) {
           (mat->table_c_rho_T[i_rho * mat->num_T + i_T] > 0)) {
         mat->c_tiny = mat->table_c_rho_T[i_rho * mat->num_T + i_T];
       }
+      if ((mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] < mat->s_tiny) &&
+          (mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] > 0)) {
+        mat->s_tiny = mat->table_log_s_rho_T[i_rho * mat->num_T + i_T];
+      }
     }
   }
 
-  // Tiny pressure to allow interpolation near non-positive values
+  // Tiny values to allow interpolation near non-positive values
+  mat->u_tiny *= 1e-3f;
   mat->P_tiny *= 1e-3f;
   mat->c_tiny *= 1e-3f;
+  mat->s_tiny *= 1e-3f;
+
+  // Convert sp. int. energies to log(sp. int. energy), same for sp. entropies
+  for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
+    for (int i_T = 0; i_T < mat->num_T; i_T++) {
+      // If not positive then set very small for the log
+      if (mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] <= 0) {
+        mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] = mat->u_tiny;
+      }
+
+      mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] =
+          logf(mat->table_log_u_rho_T[i_rho * mat->num_T + i_T]);
+
+      if (mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] <= 0) {
+        mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] = mat->s_tiny;
+      }
+
+      mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] =
+          logf(mat->table_log_s_rho_T[i_rho * mat->num_T + i_T]);
+    }
+  }
 }
 
 // Convert to internal units
@@ -255,7 +275,7 @@ INLINE static void convert_units_SESAME(struct SESAME_params *mat,
              units_cgs_conversion_factor(us, UNIT_CONV_DENSITY));
   }
 
-  // Sp. Int. Energies (log), pressures, and sound speeds
+  // Sp. int. energies (log), pressures, sound speeds, and sp. entropies
   for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
     for (int i_T = 0; i_T < mat->num_T; i_T++) {
       mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] += logf(
@@ -267,26 +287,118 @@ INLINE static void convert_units_SESAME(struct SESAME_params *mat,
       mat->table_c_rho_T[i_rho * mat->num_T + i_T] *=
           1e3f * units_cgs_conversion_factor(&si, UNIT_CONV_SPEED) /
           units_cgs_conversion_factor(us, UNIT_CONV_SPEED);
-      mat->table_s_rho_T[i_rho * mat->num_T + i_T] *=
-          units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS) /
-          units_cgs_conversion_factor(us, UNIT_CONV_ENTROPY);
+      mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] +=
+          logf(units_cgs_conversion_factor(
+                   &si, UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS) /
+               units_cgs_conversion_factor(
+                   us, UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS));
     }
   }
 
-  // Tiny pressure and sound speed
+  // Tiny values
+  mat->u_tiny *=
+      units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS) /
+      units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
   mat->P_tiny *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE) /
                  units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
   mat->c_tiny *= 1e3f * units_cgs_conversion_factor(&si, UNIT_CONV_SPEED) /
                  units_cgs_conversion_factor(us, UNIT_CONV_SPEED);
+  mat->s_tiny *=
+      units_cgs_conversion_factor(&si,
+                                  UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS) /
+      units_cgs_conversion_factor(us, UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS);
 }
 
 // gas_internal_energy_from_entropy
 INLINE static float SESAME_internal_energy_from_entropy(
     float density, float entropy, const struct SESAME_params *mat) {
 
-  error("This EOS function is not yet implemented!");
+  float u, log_u_1, log_u_2, log_u_3, log_u_4;
 
-  return 0.f;
+  if (entropy <= 0.f) {
+    return 0.f;
+  }
+
+  int idx_rho, idx_s_1, idx_s_2;
+  float intp_rho, intp_s_1, intp_s_2;
+  const float log_rho = logf(density);
+  const float log_s = logf(entropy);
+
+  // 2D interpolation (bilinear with log(rho), log(s)) to find u(rho, s))
+  // Density index
+  idx_rho =
+      find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
+
+  // Sp. entropy at this and the next density (in relevant slice of s array)
+  idx_s_1 = find_value_in_monot_incr_array(
+      log_s, mat->table_log_s_rho_T + idx_rho * mat->num_T, mat->num_T);
+  idx_s_2 = find_value_in_monot_incr_array(
+      log_s, mat->table_log_s_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
+
+  // If outside the table then extrapolate from the edge and edge-but-one values
+  if (idx_rho <= -1) {
+    idx_rho = 0;
+  } else if (idx_rho >= mat->num_rho) {
+    idx_rho = mat->num_rho - 2;
+  }
+  if (idx_s_1 <= -1) {
+    idx_s_1 = 0;
+  } else if (idx_s_1 >= mat->num_T) {
+    idx_s_1 = mat->num_T - 2;
+  }
+  if (idx_s_2 <= -1) {
+    idx_s_2 = 0;
+  } else if (idx_s_2 >= mat->num_T) {
+    idx_s_2 = mat->num_T - 2;
+  }
+
+  // Check for duplicates in SESAME tables before interpolation
+  if (mat->table_log_rho[idx_rho + 1] != mat->table_log_rho[idx_rho]) {
+    intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
+               (mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]);
+  } else {
+    intp_rho = 1.f;
+  }
+  if (mat->table_log_s_rho_T[idx_rho * mat->num_T + (idx_s_1 + 1)] !=
+      mat->table_log_s_rho_T[idx_rho * mat->num_T + idx_s_1]) {
+    intp_s_1 =
+        (log_s - mat->table_log_s_rho_T[idx_rho * mat->num_T + idx_s_1]) /
+        (mat->table_log_s_rho_T[idx_rho * mat->num_T + (idx_s_1 + 1)] -
+         mat->table_log_s_rho_T[idx_rho * mat->num_T + idx_s_1]);
+  } else {
+    intp_s_1 = 1.f;
+  }
+  if (mat->table_log_s_rho_T[(idx_rho + 1) * mat->num_T + (idx_s_2 + 1)] !=
+      mat->table_log_s_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2]) {
+    intp_s_2 =
+        (log_s - mat->table_log_s_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2]) /
+        (mat->table_log_s_rho_T[(idx_rho + 1) * mat->num_T + (idx_s_2 + 1)] -
+         mat->table_log_s_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2]);
+  } else {
+    intp_s_2 = 1.f;
+  }
+
+  // Table values
+  log_u_1 = mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_s_1];
+  log_u_2 = mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_s_1 + 1];
+  log_u_3 = mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2];
+  log_u_4 = mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2 + 1];
+
+  // If below the minimum s at this rho then just use the lowest table values
+  if ((idx_rho > 0.f) && ((intp_s_1 < 0.f) || (intp_s_2 < 0.f) ||
+                          (log_u_1 > log_u_2) || (log_u_3 > log_u_4))) {
+    intp_s_1 = 0;
+    intp_s_2 = 0;
+  }
+
+  // Interpolate with the log values
+  u = (1.f - intp_rho) * ((1.f - intp_s_1) * log_u_1 + intp_s_1 * log_u_2) +
+      intp_rho * ((1.f - intp_s_2) * log_u_3 + intp_s_2 * log_u_4);
+
+  // Convert back from log
+  u = expf(u);
+
+  return u;
 }
 
 // gas_pressure_from_entropy
@@ -338,7 +450,7 @@ INLINE static float SESAME_pressure_from_internal_energy(
   const float log_rho = logf(density);
   const float log_u = logf(u);
 
-  // 2D interpolation (bilinear with log(rho), log(u)) to find P(rho, u)
+  // 2D interpolation (bilinear with log(rho), log(u)) to find P(rho, u))
   // Density index
   idx_rho =
       find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
@@ -463,7 +575,7 @@ INLINE static float SESAME_soundspeed_from_internal_energy(
   const float log_rho = logf(density);
   const float log_u = logf(u);
 
-  // 2D interpolation (bilinear with log(rho), log(u)) to find c(rho, u)
+  // 2D interpolation (bilinear with log(rho), log(u)) to find c(rho, u))
   // Density index
   idx_rho =
       find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 0ee99bf33525317128df70869ae1c5b32d87bc3b..8dab9acba8e1ef804b11e0ef9bb56c61eb9c1928 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -25,7 +25,7 @@
  * equations)
  *
  * The thermal variable is the internal energy (u). Simple constant
- * viscosity term without switches is implemented. No thermal conduction
+ * viscosity term with the Balsara (1995) switch. No thermal conduction
  * term is implemented.
  *
  * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
diff --git a/src/hydro/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h
index 514fcad2cb1ee646d0f69dabf9c79fbaf1e8427a..f66578cf85cf0c4602005905680e82fb618f2e29 100644
--- a/src/hydro/Minimal/hydro_debug.h
+++ b/src/hydro/Minimal/hydro_debug.h
@@ -24,7 +24,7 @@
  * @brief Minimal conservative implementation of SPH (Debugging routines)
  *
  * The thermal variable is the internal energy (u). Simple constant
- * viscosity term without switches is implemented. No thermal conduction
+ * viscosity term with the Balsara (1995) switch. No thermal conduction
  * term is implemented.
  *
  * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 0c7e136d3042980f8e34698d81a3416805bb3d92..768b0358ab3457813fe503f77d69bb7828381af4 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -24,7 +24,7 @@
  * @brief Minimal conservative implementation of SPH (Neighbour loop equations)
  *
  * The thermal variable is the internal energy (u). Simple constant
- * viscosity term without switches is implemented. No thermal conduction
+ * viscosity term with the Balsara (1995) switch. No thermal conduction
  * term is implemented.
  *
  * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index dbd4d6381fadf79544d0f3536cd322d7269bc263..f5d6d51ac5d57c0789e4c140494d492e8a17bdde 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/src/hydro/Minimal/hydro_io.h
@@ -24,7 +24,7 @@
  * @brief Minimal conservative implementation of SPH (i/o routines)
  *
  * The thermal variable is the internal energy (u). Simple constant
- * viscosity term without switches is implemented. No thermal conduction
+ * viscosity term with the Balsara (1995) switch. No thermal conduction
  * term is implemented.
  *
  * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
@@ -225,8 +225,9 @@ INLINE static 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",
-                       "Minimal treatment as in Monaghan (1992)");
+  io_write_attribute_s(
+      h_grpsph, "Viscosity Model",
+      "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
 }
 
 /**
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index 158c42415e320d4f7893f8b26dd7c0296667eba0..6db84cdc0179b42d154fd96e498823f6e3f41fac 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -24,7 +24,7 @@
  * @brief Minimal conservative implementation of SPH (Particle definition)
  *
  * The thermal variable is the internal energy (u). Simple constant
- * viscosity term without switches is implemented. No thermal conduction
+ * viscosity term with the Balsara (1995) switch. No thermal conduction
  * term is implemented.
  *
  * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
diff --git a/src/hydro/Planetary/hydro.h b/src/hydro/Planetary/hydro.h
index 6b0357080527bf54065212a15516a4515df5d459..85c02aafc27a0656c05acda78dfcb75757e17b1b 100644
--- a/src/hydro/Planetary/hydro.h
+++ b/src/hydro/Planetary/hydro.h
@@ -392,10 +392,15 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
   /* Now recompute the extra quantities */
 
   /* Compute the sound speed */
+  const float pressure =
+      gas_pressure_from_internal_energy(p->rho, p->u, p->mat_id);
   const float soundspeed = hydro_get_comoving_soundspeed(p);
 
   /* Update variables. */
+  p->force.pressure = pressure;
   p->force.soundspeed = soundspeed;
+
+  p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
 }
 
 /**
@@ -472,6 +477,10 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
   p->density.wcount_dh = 0.f;
   p->rho = 0.f;
   p->density.rho_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;
 }
 
 /**
@@ -507,6 +516,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   p->density.rho_dh *= h_inv_dim_plus_one;
   p->density.wcount *= h_inv_dim;
   p->density.wcount_dh *= h_inv_dim_plus_one;
+
+  const float rho_inv = 1.f / p->rho;
+  const float a_inv2 = cosmo->a2_inv;
+
+  /* Finish calculation of the (physical) velocity curl components */
+  p->density.rot_v[0] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+  p->density.rot_v[1] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+  p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+
+  /* Finish calculation of the (physical) velocity divergence */
+  p->density.div_v *= h_inv_dim_plus_one * a_inv2 * rho_inv;
 }
 
 /**
@@ -534,6 +554,10 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
   p->density.wcount = kernel_root * h_inv_dim;
   p->density.rho_dh = 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;
 }
 
 /**
@@ -558,15 +582,22 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     const struct cosmology *cosmo, const struct hydro_props *hydro_props,
     const float dt_alpha) {
 
-  const float fac_mu = cosmo->a_factor_mu;
+  const float fac_Balsara_eps = cosmo->a_factor_Balsara_eps;
 
   /* 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 norm of div v including the Hubble flow term */
+  const float div_physical_v = p->density.div_v + hydro_dimension * cosmo->H;
+  const float abs_div_physical_v = fabsf(div_physical_v);
+
+#ifdef PLANETARY_FIXED_ENTROPY
+  /* Override the internal energy to satisfy the fixed entropy */
+  p->u = gas_internal_energy_from_entropy(p->rho, p->s_fixed, p->mat_id);
+  xp->u_full = p->u;
+#endif
 
   /* Compute the pressure */
   const float pressure =
@@ -576,23 +607,30 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float soundspeed =
       gas_soundspeed_from_internal_energy(p->rho, p->u, p->mat_id);
 
-  /* Compute the "grad h" term */
-  const float rho_inv = 1.f / p->rho;
-  float rho_dh = p->density.rho_dh;
+  /* Compute the "grad h" term  - Note here that we have \tilde{x}
+   * as 1 as we use the local number density to find neighbours. This
+   * introduces a j-component that is considered in the force loop,
+   * meaning that this cached grad_h_term gives:
+   *
+   * f_ij = 1.f - grad_h_term_i / m_j */
+  const float common_factor = p->h * hydro_dimension_inv / p->density.wcount;
+  float grad_h_term = common_factor * p->density.rho_dh /
+                      (1.f + common_factor * p->density.wcount_dh);
+
   /* Ignore changing-kernel effects when h ~= h_max */
   if (p->h > 0.9999f * hydro_props->h_max) {
-    rho_dh = 0.f;
+    grad_h_term = 0.f;
   }
-  const float grad_h_term =
-      1.f / (1.f + hydro_dimension_inv * p->h * rho_dh * rho_inv);
 
   /* Compute the Balsara switch */
 #ifdef PLANETARY_SPH_NO_BALSARA
   const float balsara = hydro_props->viscosity.alpha;
 #else
-  const float balsara =
-      hydro_props->viscosity.alpha * abs_div_v /
-      (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
+  /* Pre-multiply in the AV factor; hydro_props are not passed to the iact
+   * functions */
+  const float balsara = hydro_props->viscosity.alpha * abs_div_physical_v /
+                        (abs_div_physical_v + curl_v +
+                         0.0001f * fac_Balsara_eps * soundspeed / p->h);
 #endif
 
   /* Update variables. */
@@ -641,12 +679,12 @@ __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 */
+  /* Re-set the internal energy */
   p->u = xp->u_full;
 
   /* Compute the pressure */
   const float pressure =
-      gas_pressure_from_internal_energy(p->rho, xp->u_full, p->mat_id);
+      gas_pressure_from_internal_energy(p->rho, p->u, p->mat_id);
 
   /* Compute the sound speed */
   const float soundspeed =
@@ -654,6 +692,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
 
   p->force.pressure = pressure;
   p->force.soundspeed = soundspeed;
+
+  p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
 }
 
 /**
@@ -683,7 +723,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
   p->u += p->u_dt * dt_therm;
 
   /* Check against absolute minimum */
-  const float min_u = hydro_props->minimal_internal_energy;
+  const float min_u =
+      hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
 
   p->u = max(p->u, min_u);
 
@@ -713,6 +754,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 
   p->force.pressure = pressure;
   p->force.soundspeed = soundspeed;
+
+  p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
 }
 
 /**
@@ -759,7 +802,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
   xp->u_full = max(xp->u_full + delta_u, 0.5f * xp->u_full);
 
   /* Check against absolute minimum */
-  const float min_u = hydro_props->minimal_internal_energy;
+  const float min_u =
+      hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
 
   if (xp->u_full < min_u) {
     xp->u_full = min_u;
diff --git a/src/hydro/Planetary/hydro_iact.h b/src/hydro/Planetary/hydro_iact.h
index 426a1672e66e94067be153682534f03c31bfd8ee..b9a22c8382da6c32608cdf4194ceb015fc0bc7e6 100644
--- a/src/hydro/Planetary/hydro_iact.h
+++ b/src/hydro/Planetary/hydro_iact.h
@@ -56,6 +56,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 
   float wi, wj, wi_dx, wj_dx;
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (pi->time_bin >= time_bin_inhibited)
+    error("Inhibited pi in interaction function!");
+  if (pj->time_bin >= time_bin_inhibited)
+    error("Inhibited pj in interaction function!");
+#endif
+
   /* Get r. */
   const float r_inv = 1.0f / sqrtf(r2);
   const float r = r2 * r_inv;
@@ -83,6 +90,33 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
   pj->density.wcount += wj;
   pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
+
+  /* Compute dv dot r */
+  float dv[3], curlvr[3];
+
+  const float faci = mj * wi_dx * r_inv;
+  const float facj = mi * wj_dx * r_inv;
+
+  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 -= faci * dvdr;
+  pj->density.div_v -= facj * 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] += faci * curlvr[0];
+  pi->density.rot_v[1] += faci * curlvr[1];
+  pi->density.rot_v[2] += faci * curlvr[2];
+
+  pj->density.rot_v[0] += facj * curlvr[0];
+  pj->density.rot_v[1] += facj * curlvr[1];
+  pj->density.rot_v[2] += facj * curlvr[2];
 }
 
 /**
@@ -103,6 +137,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
 
   float wi, wi_dx;
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (pi->time_bin >= time_bin_inhibited)
+    error("Inhibited pi in interaction function!");
+  if (pj->time_bin >= time_bin_inhibited)
+    error("Inhibited pj in interaction function!");
+#endif
+
   /* Get the masses. */
   const float mj = pj->mass;
 
@@ -118,6 +159,27 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
   pi->density.wcount += wi;
   pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+  /* Compute dv dot r */
+  float dv[3], curlvr[3];
+
+  const float faci = mj * wi_dx * r_inv;
+
+  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 -= faci * 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] += faci * curlvr[0];
+  pi->density.rot_v[1] += faci * curlvr[1];
+  pi->density.rot_v[2] += faci * curlvr[2];
 }
 
 /**
@@ -136,6 +198,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
     float r2, const float *dx, float hi, float hj, struct part *restrict pi,
     struct part *restrict pj, float a, float H) {
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (pi->time_bin >= time_bin_inhibited)
+    error("Inhibited pi in interaction function!");
+  if (pj->time_bin >= time_bin_inhibited)
+    error("Inhibited pj in interaction function!");
+#endif
+
   /* Cosmological factors entering the EoMs */
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
@@ -168,9 +237,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   kernel_deval(xj, &wj, &wj_dx);
   const float wj_dr = hjd_inv * wj_dx;
 
+  /* Variable smoothing length term */
+  const float f_ij = 1.f - pi->force.f / mj;
+  const float f_ji = 1.f - pj->force.f / mi;
+
   /* Compute gradient terms */
-  const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
-  const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
+  const float P_over_rho2_i = pressurei / (rhoi * rhoi) * f_ij;
+  const float P_over_rho2_j = pressurej / (rhoj * rhoj) * f_ji;
 
   /* Compute dv dot r. */
   const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
@@ -195,7 +268,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
   /* Convolve with the kernel */
-  const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+  const float visc_acc_term =
+      0.5f * visc * (wi_dr * f_ij + wj_dr * f_ji) * r_inv;
 
   /* SPH acceleration term */
   const float sph_acc_term =
@@ -253,6 +327,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
     float r2, const float *dx, float hi, float hj, struct part *restrict pi,
     const struct part *restrict pj, float a, float H) {
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (pi->time_bin >= time_bin_inhibited)
+    error("Inhibited pi in interaction function!");
+  if (pj->time_bin >= time_bin_inhibited)
+    error("Inhibited pj in interaction function!");
+#endif
+
   /* Cosmological factors entering the EoMs */
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
@@ -262,7 +343,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float r = r2 * r_inv;
 
   /* Recover some data */
-  // const float mi = pi->mass;
+  const float mi = pi->mass;
   const float mj = pj->mass;
   const float rhoi = pi->rho;
   const float rhoj = pj->rho;
@@ -285,9 +366,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   kernel_deval(xj, &wj, &wj_dx);
   const float wj_dr = hjd_inv * wj_dx;
 
+  /* Variable smoothing length term */
+  const float f_ij = 1.f - pi->force.f / mj;
+  const float f_ji = 1.f - pj->force.f / mi;
+
   /* Compute gradient terms */
-  const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
-  const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
+  const float P_over_rho2_i = pressurei / (rhoi * rhoi) * f_ij;
+  const float P_over_rho2_j = pressurej / (rhoj * rhoj) * f_ji;
 
   /* Compute dv dot r. */
   const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
@@ -314,7 +399,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
   /* Convolve with the kernel */
-  const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+  const float visc_acc_term =
+      0.5f * visc * (wi_dr * f_ij + wj_dr * f_ji) * r_inv;
 
   /* SPH acceleration term */
   const float sph_acc_term =
diff --git a/src/hydro/Planetary/hydro_io.h b/src/hydro/Planetary/hydro_io.h
index 9872441e2d5a641069e589154ed4b52600db3f45..eabf77ebb83ed974208e559e0910e387a94ab555 100644
--- a/src/hydro/Planetary/hydro_io.h
+++ b/src/hydro/Planetary/hydro_io.h
@@ -51,7 +51,11 @@ INLINE static void hydro_read_particles(struct part* parts,
                                         struct io_props* list,
                                         int* num_fields) {
 
+#ifdef PLANETARY_FIXED_ENTROPY
+  *num_fields = 10;
+#else
   *num_fields = 9;
+#endif
 
   /* List what we want to read */
   list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
@@ -72,6 +76,11 @@ INLINE static void hydro_read_particles(struct part* parts,
                                 UNIT_CONV_DENSITY, parts, rho);
   list[8] = io_make_input_field("MaterialIDs", INT, 1, COMPULSORY,
                                 UNIT_CONV_NO_UNITS, parts, mat_id);
+#ifdef PLANETARY_FIXED_ENTROPY
+  list[9] = io_make_input_field("Entropies", FLOAT, 1, COMPULSORY,
+                                UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS, parts,
+                                s_fixed);
+#endif
 }
 
 INLINE static void convert_S(const struct engine* e, const struct part* p,
diff --git a/src/hydro/Planetary/hydro_part.h b/src/hydro/Planetary/hydro_part.h
index 482f44288c9218fb55bf6d199978b423b84c6caa..bcd905acc4a60b1d7f05bd6860b1597dc9c17df7 100644
--- a/src/hydro/Planetary/hydro_part.h
+++ b/src/hydro/Planetary/hydro_part.h
@@ -208,6 +208,11 @@ struct part {
 
 #endif
 
+#ifdef PLANETARY_FIXED_ENTROPY
+  /* Fixed specific entropy */
+  float s_fixed;
+#endif
+
 } SWIFT_STRUCT_ALIGN;
 
 #endif /* SWIFT_PLANETARY_HYDRO_PART_H */
diff --git a/src/units.c b/src/units.c
index 00c53d89953b38d09de876aea6cccf167b79f8cb..074b553c3da6f0c5209dcc2a98a528c9f70ec838 100644
--- a/src/units.c
+++ b/src/units.c
@@ -329,6 +329,12 @@ void units_get_base_unit_exponents_array(float baseUnitsExp[5],
       baseUnitsExp[UNIT_TIME] = -2.f;
       break;
 
+    case UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS:
+      baseUnitsExp[UNIT_LENGTH] = 2.f;
+      baseUnitsExp[UNIT_TIME] = -2.f;
+      baseUnitsExp[UNIT_TEMPERATURE] = -1.f;
+      break;
+
     case UNIT_CONV_POWER:
       baseUnitsExp[UNIT_MASS] = 1.f;
       baseUnitsExp[UNIT_LENGTH] = 2.f;
diff --git a/src/units.h b/src/units.h
index 87ab963e0ac640521501475121fa77ad1faaf2bc..933ab6a66c69511d43a901739a407ff7e14d3a43 100644
--- a/src/units.h
+++ b/src/units.h
@@ -84,6 +84,7 @@ enum unit_conversion_factor {
   UNIT_CONV_ENERGY_PER_UNIT_MASS,
   UNIT_CONV_ENTROPY,
   UNIT_CONV_ENTROPY_PER_UNIT_MASS,
+  UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS,
   UNIT_CONV_POWER,
   UNIT_CONV_PRESSURE,
   UNIT_CONV_FREQUENCY,