diff --git a/src/hydro/SPHENIX/hydro.h b/src/hydro/SPHENIX/hydro.h
index a344d982cd0a01586e23de47d2006d904ffac196..773f8e39a1819017e1433eb9cd7bca549aaf849d 100644
--- a/src/hydro/SPHENIX/hydro.h
+++ b/src/hydro/SPHENIX/hydro.h
@@ -407,7 +407,7 @@ hydro_set_drifted_physical_internal_energy(
p->force.soundspeed = soundspeed;
p->force.pressure = pressure_including_floor;
- p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
+ p->viscosity.v_sig = max(p->viscosity.v_sig, soundspeed);
}
/**
@@ -499,13 +499,11 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* @brief beta The non-linear viscosity constant.
*/
__attribute__((always_inline)) INLINE static float hydro_signal_velocity(
- const float dx[3], const struct part *restrict pi,
- const struct part *restrict pj, const float mu_ij, const float beta) {
+ const struct part *restrict pi, const float alphai, const float mu_ij, const float beta) {
const float ci = pi->force.soundspeed;
- const float cj = pj->force.soundspeed;
- return ci + cj - beta * mu_ij;
+ return alphai * ci + beta * fabsf(mu_ij);
}
/**
@@ -720,7 +718,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_gradient(
p->viscosity.div_v = 0.f;
- p->viscosity.v_sig = 2.f * p->force.soundspeed;
+ p->viscosity.v_sig = p->force.soundspeed;
p->force.alpha_visc_max_ngb = p->viscosity.alpha;
}
@@ -1015,7 +1013,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
p->force.soundspeed = soundspeed;
/* Update the signal velocity, if we need to. */
- p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
+ p->viscosity.v_sig = max(p->viscosity.v_sig, soundspeed);
}
/**
@@ -1090,7 +1088,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->force.soundspeed = soundspeed;
/* Update signal velocity if we need to */
- p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
+ p->viscosity.v_sig = max(p->viscosity.v_sig, soundspeed);
}
/**
diff --git a/src/hydro/SPHENIX/hydro_iact.h b/src/hydro/SPHENIX/hydro_iact.h
index 197bf110fcbbf232b45152ce971d03da7c71ac64..3dce5424dc770f982b8a5aac9308e449d539c7c8 100644
--- a/src/hydro/SPHENIX/hydro_iact.h
+++ b/src/hydro/SPHENIX/hydro_iact.h
@@ -211,12 +211,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Signal velocity */
- const float new_v_sig =
- signal_velocity(pi, pj, mu_ij, const_viscosity_beta);
+ const float new_v_sig_i =
+ signal_velocity(pi, 1.0f, mu_ij, const_viscosity_beta);
+ const float new_v_sig_j =
+ signal_velocity(pj, 1.0f, mu_ij, const_viscosity_beta);
/* Update if we need to */
- pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig);
- pj->viscosity.v_sig = max(pj->viscosity.v_sig, new_v_sig);
+ pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig_i);
+ pj->viscosity.v_sig = max(pj->viscosity.v_sig, new_v_sig_j);
/* Now we need to compute the div terms */
const float faci = mj * f_ij * wi_dx * r_inv;
@@ -321,11 +323,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Signal velocity */
- const float new_v_sig =
- signal_velocity(pi, pj, mu_ij, const_viscosity_beta);
+ const float new_v_sig_i =
+ signal_velocity(pi, 1.0f, mu_ij, const_viscosity_beta);
/* Update if we need to */
- pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig);
+ pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig_i);
/* Now we need to compute the div terms */
const float faci = mj * f_ij * wi_dx * r_inv;
@@ -420,9 +422,22 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float omega_ij = min(dvdr_Hubble, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+ /* Update timestep if needws be */
+ const float v_sig_dt_i =
+ signal_velocity(pi, 1.0f, mu_ij, const_viscosity_beta);
+ const float v_sig_dt_j =
+ signal_velocity(pj, 1.0f, mu_ij, const_viscosity_beta);
+
+ pi->viscosity.v_sig = max(pi->viscosity.v_sig, v_sig_dt_i);
+ pj->viscosity.v_sig = max(pj->viscosity.v_sig, v_sig_dt_j);
+
/* Compute sound speeds and signal velocity */
- const float v_sig =
- signal_velocity(pi, pj, mu_ij, const_viscosity_beta);
+ const float alphai = pi->viscosity.alpha;
+ const float alphaj = pj->viscosity.alpha;
+ const float v_sig_i =
+ signal_velocity(pi, alphai, mu_ij, const_viscosity_beta);
+ const float v_sig_j =
+ signal_velocity(pj, alphaj, mu_ij, const_viscosity_beta);
/* Variable smoothing length term */
const float f_ij = 1.f - pi->force.f / mj;
@@ -432,16 +447,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float balsara_i = pi->force.balsara;
const float balsara_j = pj->force.balsara;
+ const float q_ij = -0.5f * rhoi * v_sig_i * mu_ij;
+ const float q_ji = -0.5f * rhoj * v_sig_j * mu_ij;
+
/* Construct the full viscosity term */
- const float rho_ij = rhoi + rhoj;
- const float alpha = pi->viscosity.alpha + pj->viscosity.alpha;
- const float visc =
- -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
-
- /* Convolve with the kernel */
- const float visc_acc_term =
- 0.5f * visc * (wi_dr * f_ij + wj_dr * f_ji) * r_inv;
-
+ float visc_acc_term = f_ij * balsara_i * q_ij * wi_dr / (rhoi * rhoi) * r_inv;
+ visc_acc_term += f_ji * balsara_j * q_ji * wj_dr / (rhoj * rhoj) * r_inv;
+
/* Compute gradient terms */
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * f_ij;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * f_ji;
@@ -478,6 +490,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
* alpha from the highest pressure particle to dominate, so that the
* diffusion limited particles always take precedence - another trick to
* allow the scheme to work with thermal feedback. */
+ const float rho_ij = rhoi + rhoj;
const float alpha_diff =
(pressurei * pi->diffusion.alpha + pressurej * pj->diffusion.alpha) /
(pressurei + pressurej);
@@ -572,10 +585,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float omega_ij = min(dvdr_Hubble, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+ /* Update timestep if needws be */
+ const float v_sig_dt_i =
+ signal_velocity(pi, 1.0f, mu_ij, const_viscosity_beta);
+
+ pi->viscosity.v_sig = max(pi->viscosity.v_sig, v_sig_dt_i);
+
/* Compute sound speeds and signal velocity */
- const float v_sig =
- signal_velocity(pi, pj, mu_ij, const_viscosity_beta);
-
+ const float alphai = pi->viscosity.alpha;
+ const float alphaj = pj->viscosity.alpha;
+ const float v_sig_i =
+ signal_velocity(pi, alphai, mu_ij, const_viscosity_beta);
+ const float v_sig_j =
+ signal_velocity(pj, alphaj, mu_ij, const_viscosity_beta);
+
/* Variable smoothing length term */
const float f_ij = 1.f - pi->force.f / mj;
const float f_ji = 1.f - pj->force.f / mi;
@@ -584,16 +607,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float balsara_i = pi->force.balsara;
const float balsara_j = pj->force.balsara;
- /* Construct the full viscosity term */
- const float rho_ij = rhoi + rhoj;
- const float alpha = pi->viscosity.alpha + pj->viscosity.alpha;
- const float visc =
- -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
-
- /* Convolve with the kernel */
- const float visc_acc_term =
- 0.5f * visc * (wi_dr * f_ij + wj_dr * f_ji) * r_inv;
+ const float q_ij = -0.5f * rhoi * v_sig_i * mu_ij;
+ const float q_ji = -0.5f * rhoj * v_sig_j * mu_ij;
+ /* Construct the full viscosity term */
+ float visc_acc_term = f_ij * balsara_i * q_ij * wi_dr / (rhoi * rhoi) * r_inv;
+ visc_acc_term += f_ji * balsara_j * q_ji * wj_dr / (rhoj * rhoj) * r_inv;
+
/* Compute gradient terms */
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * f_ij;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * f_ji;
@@ -625,6 +645,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
* alpha from the highest pressure particle to dominate, so that the
* diffusion limited particles always take precedence - another trick to
* allow the scheme to work with thermal feedback. */
+ const float rho_ij = rhoi + rhoj;
const float alpha_diff =
(pressurei * pi->diffusion.alpha + pressurej * pj->diffusion.alpha) /
(pressurei + pressurej);
diff --git a/src/hydro/SPHENIX/hydro_parameters.h b/src/hydro/SPHENIX/hydro_parameters.h
index 277be9562dd079d160ca1a9d5ab7208551c1e276..103b9d1628a45a57fd35ac6776f0fd50a294119b 100644
--- a/src/hydro/SPHENIX/hydro_parameters.h
+++ b/src/hydro/SPHENIX/hydro_parameters.h
@@ -50,7 +50,7 @@
/*! Cosmology default beta=3.0.
* Alpha can be set in the parameter file.
* Beta is defined as in e.g. Price (2010) Eqn (103) */
-#define const_viscosity_beta 3.0f
+#define const_viscosity_beta 2.0f
/*! The viscosity that the particles are reset to after being hit by a
* feedback event. This should be set to the same value as the
@@ -68,7 +68,7 @@
#define hydro_props_default_viscosity_alpha_min 0.0f
/*! Maximal value for the viscosity alpha in variable schemes. */
-#define hydro_props_default_viscosity_alpha_max 2.0f
+#define hydro_props_default_viscosity_alpha_max 1.0f
/*! Decay length for the viscosity scheme. This is scheme dependent. */
#define hydro_props_default_viscosity_length 0.05f
diff --git a/src/mhd/DirectInduction/mhd.h b/src/mhd/DirectInduction/mhd.h
index 58f5a48398e63896c00ca3ee37241d1f79b1b959..2f956c9c5a6e9eca297a850569292fd85b407e97 100644
--- a/src/mhd/DirectInduction/mhd.h
+++ b/src/mhd/DirectInduction/mhd.h
@@ -247,12 +247,11 @@ mhd_get_fast_magnetosonic_wave_phase_velocity(const float dx[3], const struct pa
* @brief beta The non-linear viscosity constant.
*/
__attribute__((always_inline)) INLINE static float mhd_signal_velocity(
- const struct part *restrict pi, const struct part *restrict pj, const float mu_ij, const float beta) {
+ const struct part *restrict pi, const float alphai, const float mu_ij, const float beta) {
const float cmsi = mhd_get_magnetosonic_speed(pi);
- const float cmsj = mhd_get_magnetosonic_speed(pj);
- const float v_sig = cmsi + cmsj + beta * fabsf(mu_ij);
+ const float v_sig = alphai * cmsi + beta * fabsf(mu_ij);
return v_sig;
}
diff --git a/src/signal_velocity.h b/src/signal_velocity.h
index 7e476b7b79dcef95a471e7c2738a887c78f0de9a..160eaf6442ff6645c035b801a0b0b6af7199ae04 100644
--- a/src/signal_velocity.h
+++ b/src/signal_velocity.h
@@ -41,9 +41,9 @@
* @brief beta The non-linear viscosity constant.
*/
__attribute__((always_inline)) INLINE static float signal_velocity(
- const struct part *restrict pi, const struct part *restrict pj, const float mu_ij, const float beta) {
+ const struct part *restrict pi, const float alphai, const float mu_ij, const float beta) {
- return mhd_signal_velocity(pi, pj, mu_ij, beta);
+ return mhd_signal_velocity(pi, alphai, mu_ij, beta);
}
#else
@@ -65,9 +65,9 @@ __attribute__((always_inline)) INLINE static float signal_velocity(
* @brief beta The non-linear viscosity constant.
*/
__attribute__((always_inline)) INLINE static float signal_velocity(
- const struct part *restrict pi,const struct part *restrict pj, const float mu_ij, const float beta) {
+ const struct part *restrict pi, const float alphai, const float mu_ij, const float beta) {
- return hydro_signal_velocity(pi, pj, mu_ij, beta);
+ return hydro_signal_velocity(pi, alphai, mu_ij, beta);
}
#endif