diff --git a/src/star_formation/GEAR/star_formation.h b/src/star_formation/GEAR/star_formation.h index cfb314be4384b21d686e2bdca4b11521ccb978b0..a1d9dda2a0736bde67167f4e11e2de4fcb8c72cc 100644 --- a/src/star_formation/GEAR/star_formation.h +++ b/src/star_formation/GEAR/star_formation.h @@ -85,11 +85,7 @@ INLINE static int star_formation_is_star_forming( / (mu * phys_const->const_proton_mass) + sigma2); /* Check the density criterion */ - if (density > density_criterion) { - return 1; - } else { - return 0; - } + return density > density_criterion; } /** @@ -139,20 +135,15 @@ INLINE static int star_formation_should_convert_to_star( const float density = hydro_get_physical_density(p, cosmo); /* Compute the probability */ - const float inv_free_fall_time = sqrtf(density * 32. * G / (3. * M_PI)); - const float prob = 1. - exp(-starform->star_formation_efficiency * inv_free_fall_time * dt_star); + const float inv_free_fall_time = sqrtf(density * 32.f * G * 0.33333333f * M_1_PI); + const float prob = 1.f - exp(-starform->star_formation_efficiency * inv_free_fall_time * dt_star); /* Roll the dice... */ const float random_number = random_unit_interval(p->id, e->ti_current, random_number_star_formation); - if (random_number > prob) { - /* No star for you */ - return 0; - } else { - /* You get a star, you get a star, everybody gets a star */ - return 1; - } + /* Can we form a star? */ + return random_number < prob; } /** @@ -241,6 +232,9 @@ __attribute__((always_inline)) INLINE static void star_formation_end_density( // TODO move into pressure floor /* To finish the turbulence estimation we devide by the density */ p->sf_data.sigma2 /= pow_dimension(p->h) * hydro_get_physical_density(p, cosmo); + + /* Add the cosmological factor */ + p->sf_data.sigma2 *= cosmo->a * cosmo->a; } /** diff --git a/src/star_formation/GEAR/star_formation_iact.h b/src/star_formation/GEAR/star_formation_iact.h index ef42d264ee6f870694c185015d1ad2bc15c72223..80de6422757e8f7d22e5691a70e5e9d21adb09b1 100644 --- a/src/star_formation/GEAR/star_formation_iact.h +++ b/src/star_formation/GEAR/star_formation_iact.h @@ -57,12 +57,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_star_formation( pi->v[2] - pj->v[2] }; - /* Square of delta v */ - float norm_v2 = dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2]; + /* Norms at power 2 */ + const float norm_v2 = dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2]; + const float norm_x2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; /* Compute the velocity dispersion */ - pi->sf_data.sigma2 += norm_v2 * wi * hydro_get_mass(pj); - pj->sf_data.sigma2 += norm_v2 * wj * hydro_get_mass(pi); + const float sigma2 = norm_v2 + H * norm_x2; + + /* Compute the velocity dispersion */ + pi->sf_data.sigma2 += sigma2 * wi * hydro_get_mass(pj); + pj->sf_data.sigma2 += sigma2 * wj * hydro_get_mass(pi); } /** @@ -94,11 +98,15 @@ runner_iact_nonsym_star_formation(float r2, const float *dx, float hi, float hj, pi->v[2] - pj->v[2] }; - /* Square of delta v */ - float norm_v2 = dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2]; + /* Norms at power 2 */ + const float norm_v2 = dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2]; + const float norm_x2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; + + /* Compute the velocity dispersion */ + const float sigma2 = norm_v2 + H * norm_x2; /* Compute the velocity dispersion */ - pi->sf_data.sigma2 += norm_v2 * wi * hydro_get_mass(pj); + pi->sf_data.sigma2 += sigma2 * wi * hydro_get_mass(pj); } #endif /* SWIFT_GEAR_STAR_FORMATION_IACT_H */