diff --git a/src/const.h b/src/const.h
index b1d8500bcaaae97c970719af01cf685b42e57ac9..bff000ecba735aa2ca237ccebee526e99c8853fa 100644
--- a/src/const.h
+++ b/src/const.h
@@ -21,10 +21,10 @@
#define SWIFT_CONST_H
/* SPH Viscosity constants. */
-/* Cosmology defaults: a=0.8, b=1.5. Planetary defaults: a=1.5, b=2.0
- * Note: beta is multiplied by alpha, as in e.g. Springel (2010) Eqn (33) */
+/* Cosmology defaults: a=0.8, b=3.0. Planetary defaults: a=1.5, b=4.0
+ * Beta is defined as in e.g. Price (2010) Eqn (103) */
#define const_viscosity_alpha 0.8f
-#define const_viscosity_beta 1.5f
+#define const_viscosity_beta 3.0f
#define const_viscosity_alpha_min \
0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index c903f5e692fd9876b99782f20582b7bc6ff1f637..cbe650d726c1e8cb6f877e2e514d5243c41ef8a3 100644
--- a/src/hydro/Default/hydro_iact.h
+++ b/src/hydro/Default/hydro_iact.h
@@ -227,7 +227,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Compute signal velocity */
v_sig = pi->force.soundspeed + pj->force.soundspeed -
- 2.f * const_viscosity_beta * omega_ij;
+ const_viscosity_beta * omega_ij;
/* Compute viscosity parameter */
alpha_ij = -0.5f * (pi->alpha + pj->alpha);
@@ -337,7 +337,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Compute signal velocity */
v_sig = pi->force.soundspeed + pj->force.soundspeed -
- 2.f * const_viscosity_beta * omega_ij;
+ const_viscosity_beta * omega_ij;
/* Compute viscosity parameter */
alpha_ij = -0.5f * (pi->alpha + pj->alpha);
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index e46ada9b75b89e91296b6fad43ba2093024c6164..aeabe93a9f6bf2dd5742fb62634f335f31e37eb8 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -493,7 +493,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Signal velocity */
- const float v_sig = ci + cj - 2.f * const_viscosity_beta * mu_ij;
+ const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
/* Now construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
@@ -616,7 +616,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Signal velocity */
- const float v_sig = ci + cj - 2.f * const_viscosity_beta * mu_ij;
+ const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
/* Now construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
@@ -741,7 +741,7 @@ runner_iact_nonsym_1_vec_force(
vec_mul(ri.v, omega_ij.v)); /* This is 0 or negative */
/* Compute signal velocity */
- v_sig.v = vec_fnma(vec_set1(2.f * const_viscosity_beta),
+ v_sig.v = vec_fnma(vec_set1(const_viscosity_beta),
mu_ij.v, vec_add(ci.v, cj.v));
/* Now construct the full viscosity term */
@@ -928,9 +928,9 @@ runner_iact_nonsym_2_vec_force(
v_fac_mu.v, vec_mul(ri_2.v, omega_ij_2.v)); /* This is 0 or negative */
/* Compute signal velocity */
- v_sig.v = vec_fnma(vec_set1(2.f * const_viscosity_beta),
+ v_sig.v = vec_fnma(vec_set1(const_viscosity_beta),
mu_ij.v, vec_add(ci.v, cj.v));
- v_sig_2.v = vec_fnma(vec_set1(2.f * const_viscosity_beta),
+ v_sig_2.v = vec_fnma(vec_set1(const_viscosity_beta),
mu_ij_2.v, vec_add(ci.v, cj_2.v));
/* Now construct the full viscosity term */
diff --git a/src/hydro/GizmoMFM/hydro_iact.h b/src/hydro/GizmoMFM/hydro_iact.h
index 6c8d001441e7450a758cc704296b397f55c783ac..38a97cbea39c1ed5c6926c911941e655e52362aa 100644
--- a/src/hydro/GizmoMFM/hydro_iact.h
+++ b/src/hydro/GizmoMFM/hydro_iact.h
@@ -267,7 +267,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
const float dvdotdx = min(dvdr, 0.0f);
/* Get the signal velocity */
- vmax -= 2.f * const_viscosity_beta * dvdotdx * r_inv;
+ vmax -= const_viscosity_beta * dvdotdx * r_inv;
/* Store the signal velocity */
pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
diff --git a/src/hydro/GizmoMFV/hydro_iact.h b/src/hydro/GizmoMFV/hydro_iact.h
index 69b27abbbd5272f6e9552774e3592b2179549585..2f73e67ea2fdcecc527de8b1af0d15731f967b9b 100644
--- a/src/hydro/GizmoMFV/hydro_iact.h
+++ b/src/hydro/GizmoMFV/hydro_iact.h
@@ -271,7 +271,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
dvdotdx = min(dvdotdx, 0.f);
/* Get the signal velocity */
- vmax -= 2.f * const_viscosity_beta * dvdotdx * r_inv;
+ vmax -= const_viscosity_beta * dvdotdx * r_inv;
/* Store the signal velocity */
pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index f0aba6f97ab0ad569f25e19297b9b28f4fc89c99..e9ad55516b09d1b28988e68f7316edbc7352be71 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -184,7 +184,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Compute sound speeds and signal velocity */
const float ci = pi->force.soundspeed;
const float cj = pj->force.soundspeed;
- const float v_sig = ci + cj - 2.f * const_viscosity_beta * mu_ij;
+ const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
/* Construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
@@ -300,7 +300,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Compute sound speeds and signal velocity */
const float ci = pi->force.soundspeed;
const float cj = pj->force.soundspeed;
- const float v_sig = ci + cj - 2.f * const_viscosity_beta * mu_ij;
+ const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
/* Construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
diff --git a/src/hydro/Planetary/hydro_iact.h b/src/hydro/Planetary/hydro_iact.h
index 906d5c86e4793dd4405f6ebec690023694f44a46..71f616311bf0c217555c77ed40fc584733471937 100644
--- a/src/hydro/Planetary/hydro_iact.h
+++ b/src/hydro/Planetary/hydro_iact.h
@@ -189,7 +189,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Compute sound speeds and signal velocity */
const float ci = pi->force.soundspeed;
const float cj = pj->force.soundspeed;
- const float v_sig = ci + cj - 2.f * const_viscosity_beta * mu_ij;
+ const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
/* Now construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
@@ -315,7 +315,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float cj = pj->force.soundspeed;
/* Signal velocity */
- const float v_sig = ci + cj - 2.f * const_viscosity_beta * mu_ij;
+ const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
/* Construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
diff --git a/src/hydro/PressureEnergy/hydro_iact.h b/src/hydro/PressureEnergy/hydro_iact.h
index 4aba8a262732b55e10f10f11d93044864df9b305..78aa4e724d5b61a125a30bd900211ea284f0db94 100644
--- a/src/hydro/PressureEnergy/hydro_iact.h
+++ b/src/hydro/PressureEnergy/hydro_iact.h
@@ -244,7 +244,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Compute sound speeds and signal velocity */
const float ci = pi->force.soundspeed;
const float cj = pj->force.soundspeed;
- const float v_sig = ci + cj - 2.f * const_viscosity_beta * mu_ij;
+ const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
/* Balsara term */
const float balsara_i = pi->force.balsara;
@@ -372,7 +372,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Compute sound speeds and signal velocity */
const float ci = pi->force.soundspeed;
const float cj = pj->force.soundspeed;
- const float v_sig = ci + cj - 2.f * const_viscosity_beta * mu_ij;
+ const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
/* Balsara term */
const float balsara_i = pi->force.balsara;
diff --git a/src/hydro/PressureEntropy/hydro_iact.h b/src/hydro/PressureEntropy/hydro_iact.h
index 64b42093c39b58fb48b182cc5aee64f67bbb4a98..f93e99137b8ffc47c495669047d25ab17062c593 100644
--- a/src/hydro/PressureEntropy/hydro_iact.h
+++ b/src/hydro/PressureEntropy/hydro_iact.h
@@ -259,7 +259,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Signal velocity */
- const float v_sig = ci + cj - 2.f * const_viscosity_beta * mu_ij;
+ const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
/* Now construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
@@ -373,7 +373,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Signal velocity */
- const float v_sig = ci + cj - 2.f * const_viscosity_beta * mu_ij;
+ const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
/* Now construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);