From 4edc55656c13baf8532f79ff2db3986d736a104f Mon Sep 17 00:00:00 2001
From: Jacob Kegerreis <jacob.kegerreis@durham.ac.uk>
Date: Sun, 7 Oct 2018 10:47:14 +0100
Subject: [PATCH] Define the viscosity beta as a constant instead of hard-coded

---
 src/const.h                      | 5 ++++-
 src/hydro/Planetary/hydro_iact.h | 4 ++--
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/const.h b/src/const.h
index 6c5b5299c0..805b1abd45 100644
--- a/src/const.h
+++ b/src/const.h
@@ -21,7 +21,10 @@
 #define SWIFT_CONST_H
 
 /* SPH Viscosity constants. */
-#define const_viscosity_alpha 0.8f
+/* Cosmology defaults: a=0.8, b=1.5a. Planetary defaults: a=1.5, b=2.0a */ 
+#define const_viscosity_alpha 1.5f
+#define const_viscosity_beta_over_alpha 2.0f
+
 #define const_viscosity_alpha_min \
   0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
 #define const_viscosity_alpha_max \
diff --git a/src/hydro/Planetary/hydro_iact.h b/src/hydro/Planetary/hydro_iact.h
index bf96034696..ec55069128 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 - 3.f * mu_ij;
+  const float v_sig = ci + cj - 2.f * const_viscosity_beta_over_alpha * 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 - 3.f * mu_ij;
+  const float v_sig = ci + cj - 2.f * const_viscosity_beta_over_alpha * mu_ij;
 
   /* Construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-- 
GitLab