Fix diffusion velocity bug

Fixes a bug in the diffusion velocity that meant that it did not scale correctly with inter particle separation or cosmology.

Very small physics change.

src/hydro/AnarchyDU/hydro_iact.h

@@ -397,12 +397,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
 
   /* Diffusion term */
-  const float v_diff =
-      max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f);
   const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
+  const float v_diff =
+      alpha_diff * sqrtf(0.5f * fabsf(pressurei - pressurej) / rho_ij) +
+      fabsf(fac_mu * r_inv * dvdr_Hubble);
   /* wi_dx + wj_dx / 2 is F_ij */
   const float diff_du_term =
-      alpha_diff * fac_mu * v_diff * (pi->u - pj->u) * (wi_dr + wj_dr) / rho_ij;
+      v_diff * (pi->u - pj->u) * (wi_dr / rhoi + wj_dr / rhoj);
 
   /* Assemble the energy equation term */
   const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;
@@ -518,12 +519,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
 
   /* Diffusion term */
-  const float v_diff =
-      max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f);
   const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
+  const float v_diff =
+      alpha_diff * sqrtf(0.5f * fabsf(pressurei - pressurej) / rho_ij) +
+      fabsf(fac_mu * r_inv * dvdr_Hubble);
   /* wi_dx + wj_dx / 2 is F_ij */
   const float diff_du_term =
-      alpha_diff * fac_mu * v_diff * (pi->u - pj->u) * (wi_dr + wj_dr) / rho_ij;
+      v_diff * (pi->u - pj->u) * (wi_dr / rhoi + wj_dr / rhoj);
 
   /* Assemble the energy equation term */
   const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;