diff --git a/examples/SantaBarbara/README b/examples/SantaBarbara/README
new file mode 100644
index 0000000000000000000000000000000000000000..973d09a438e8ff39d7c4a91c5a08d4dfc4a51a60
--- /dev/null
+++ b/examples/SantaBarbara/README
@@ -0,0 +1,16 @@
+Initital conditions for the Santa-Barbara cluster comparison project.
+These have been regenerated from the orinigal Frenk et al. 1999 paper.
+
+The cosmology is Omega_m = 1, Omega_b = 0.1, h = 0.5 and sigma_8 = 0.9.
+
+The ICs are 256^3 particles in a 64^3 Mpc^3 volume. This is about 10x
+higher resolution than in the original paper. The ICs have been
+created for Gadget and the positions and box size are hence expressed
+in h-full units (e.g. box size of 32 / h Mpc). Similarly, the peculiar
+velocitites contain an extra sqrt(a) factor.
+
+We will use SWIFT to cancel the h- and a-factors from the ICs. Gas
+particles will be generated at startup.
+
+MD5 check-sum of the ICS:
+ba9ab4f00a70d39fa601a4a59984b343 SantaBarbara.hdf5
diff --git a/examples/SantaBarbara/getIC.sh b/examples/SantaBarbara/getIC.sh
new file mode 100755
index 0000000000000000000000000000000000000000..a3073631ceedea47c8ac218a5e62529efee6fc56
--- /dev/null
+++ b/examples/SantaBarbara/getIC.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/SantaBarbara.hdf5
diff --git a/examples/SantaBarbara/santa_barbara.yml b/examples/SantaBarbara/santa_barbara.yml
new file mode 100644
index 0000000000000000000000000000000000000000..b01321c521ec28b78c97cd082f0fd520d7024c3d
--- /dev/null
+++ b/examples/SantaBarbara/santa_barbara.yml
@@ -0,0 +1,59 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.98848e43 # 10^10 Msun
+ UnitLength_in_cgs: 3.08567758e24 # 1 Mpc
+ UnitVelocity_in_cgs: 1e5 # 1 km/s
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Cosmological parameters
+Cosmology:
+ h: 0.5
+ a_begin: 0.047619048 # z_ini = 20
+ a_end: 1.0 # z_end = 0
+ Omega_m: 1.0
+ Omega_lambda: 0.0
+ Omega_b: 0.1
+
+# Parameters governing the time integration
+TimeIntegration:
+ dt_max: 0.01
+ dt_min: 1e-10
+
+Scheduler:
+ max_top_level_cells: 16
+ cell_split_size: 100
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: santabarbara
+ scale_factor_first: 0.05
+ delta_time: 1.02
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1.02
+ scale_factor_first: 0.05
+
+# Parameters for the self-gravity scheme
+Gravity:
+ eta: 0.025
+ theta: 0.5
+ comoving_softening: 0.01 # 10 kpc = 1/25 mean inter-particle separation
+ max_physical_softening: 0.00263 # 10 ckpc = 2.63 pkpc at z=2.8 (EAGLE-like evolution of softening).
+ mesh_side_length: 128
+
+# Parameters of the hydro scheme
+SPH:
+ resolution_eta: 1.2348 # "48 Ngb" with the cubic spline kernel
+ CFL_condition: 0.1
+ initial_temperature: 1200. # (1 + z_ini)^2 * 2.72K
+ minimal_temperature: 100.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./SantaBarbara.hdf5
+ cleanup_h_factors: 1 # ICs were generated for Gadget, we need to get rid of h-factors
+ cleanup_velocity_factors: 1 # ICs were generated for Gadget, we need to get rid of sqrt(a) factors in the velocity
+ generate_gas_in_ics: 1 # Generate gas particles from the DM-only ICs
+ cleanup_smoothing_lengths: 1 # Since we generate gas, make use of the (expensive) cleaning-up procedure.
\ No newline at end of file
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 4511b2d655b0e7b3293633a466c76757a0237874..b66ada21d60555eb2ace9264a75d5cdf3dfec7d5 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -515,8 +515,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
p->force.h_dt *= p->h * hydro_dimension_inv;
- p->entropy_dt = 0.5f * cosmo->a2_inv *
- gas_entropy_from_internal_energy(p->rho, p->entropy_dt);
+ p->entropy_dt = 0.5f * gas_entropy_from_internal_energy(p->rho, p->entropy_dt);
}
/**
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index b2af8909bed1780586a5130370222c9b8157d724..55f7fb2d0ade2cdd053fbcfaf3b15f84f3f5b6d2 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -479,14 +479,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Add Hubble flow */
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
/* Balsara term */
const float balsara_i = pi->force.balsara;
const float balsara_j = pj->force.balsara;
/* Are the particles moving towards each others ? */
- const float omega_ij = (dvdr < 0.f) ? dvdr : 0.f;
+ const float omega_ij = min(dvdr_Hubble, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Signal velocity */
@@ -599,14 +602,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Add Hubble flow */
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
/* Balsara term */
const float balsara_i = pi->force.balsara;
const float balsara_j = pj->force.balsara;
/* Are the particles moving towards each others ? */
- const float omega_ij = (dvdr < 0.f) ? dvdr : 0.f;
+ const float omega_ij = min(dvdr_Hubble, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Signal velocity */
@@ -671,7 +677,7 @@ runner_iact_nonsym_1_vec_force(
vector dvx, dvy, dvz;
vector xi, xj;
vector hid_inv, hjd_inv;
- vector wi_dx, wj_dx, wi_dr, wj_dr, dvdr;
+ vector wi_dx, wj_dx, wi_dr, wj_dr, dvdr, dvdr_Hubble;
vector piax, piay, piaz;
vector pih_dt;
vector v_sig;
@@ -726,11 +732,14 @@ runner_iact_nonsym_1_vec_force(
dvdr.v =
vec_fma(dvx.v, dx->v,
vec_fma(dvy.v, dy->v,
- vec_fma(dvz.v, dz->v, vec_mul(v_a2_Hubble.v, r2->v))));
+ vec_mul(dvz.v, dz->v)));
+ /* Add Hubble flow */
+ dvdr_Hubble.v = vec_add(dvdr.v, vec_mul(v_a2_Hubble.v, r2->v));
+
/* Compute the relative velocity. (This is 0 if the particles move away from
* each other and negative otherwise) */
- omega_ij.v = vec_fmin(dvdr.v, vec_setzero());
+ omega_ij.v = vec_fmin(dvdr_Hubble.v, vec_setzero());
mu_ij.v = vec_mul(v_fac_mu.v,
vec_mul(ri.v, omega_ij.v)); /* This is 0 or negative */
@@ -806,7 +815,7 @@ runner_iact_nonsym_2_vec_force(
vector dvx, dvy, dvz;
vector ui, uj;
vector hid_inv, hjd_inv;
- vector wi_dx, wj_dx, wi_dr, wj_dr, dvdr;
+ vector wi_dx, wj_dx, wi_dr, wj_dr, dvdr, dvdr_Hubble;
vector piax, piay, piaz;
vector pih_dt;
vector v_sig;
@@ -817,7 +826,7 @@ runner_iact_nonsym_2_vec_force(
vector dvx_2, dvy_2, dvz_2;
vector ui_2, uj_2;
vector hjd_inv_2;
- vector wi_dx_2, wj_dx_2, wi_dr_2, wj_dr_2, dvdr_2;
+ vector wi_dx_2, wj_dx_2, wi_dr_2, wj_dr_2, dvdr_2, dvdr_Hubble_2;
vector piax_2, piay_2, piaz_2;
vector pih_dt_2;
vector v_sig_2;
@@ -905,16 +914,19 @@ runner_iact_nonsym_2_vec_force(
/* Compute dv dot r. */
dvdr.v = vec_fma(
dvx.v, dx.v,
- vec_fma(dvy.v, dy.v, vec_fma(dvz.v, dz.v, vec_mul(v_a2_Hubble.v, r2.v))));
+ vec_fma(dvy.v, dy.v, vec_mul(dvz.v, dz.v)));
dvdr_2.v = vec_fma(
dvx_2.v, dx_2.v,
- vec_fma(dvy_2.v, dy_2.v,
- vec_fma(dvz_2.v, dz_2.v, vec_mul(v_a2_Hubble.v, r2_2.v))));
+ vec_fma(dvy_2.v, dy_2.v, vec_mul(dvz_2.v, dz_2.v)));
+ /* Add the Hubble flow */
+ dvdr_Hubble.v = vec_add(dvdr.v, vec_mul(v_a2_Hubble.v, r2.v));
+ dvdr_Hubble_2.v = vec_add(dvdr_2.v, vec_mul(v_a2_Hubble.v, r2_2.v));
+
/* Compute the relative velocity. (This is 0 if the particles move away from
* each other and negative otherwise) */
- omega_ij.v = vec_fmin(dvdr.v, vec_setzero());
- omega_ij_2.v = vec_fmin(dvdr_2.v, vec_setzero());
+ omega_ij.v = vec_fmin(dvdr_Hubble.v, vec_setzero());
+ omega_ij_2.v = vec_fmin(dvdr_Hubble_2.v, vec_setzero());
mu_ij.v = vec_mul(v_fac_mu.v,
vec_mul(ri.v, omega_ij.v)); /* This is 0 or negative */
mu_ij_2.v = vec_mul(