Merge branch 'santa_barbara' into 'master'

Santa barbara cluster and correct cosmological terms in Gadget-SPH

See merge request !603

examples/SantaBarbara/README

+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

examples/SantaBarbara/getIC.sh

+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/SantaBarbara.hdf5

examples/SantaBarbara/santa_barbara.yml

+# 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

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);
 }
 
 /**

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(