From b885c0c4479de8d0b1fb451f855b95ff62365a25 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Mon, 29 Jul 2019 13:14:50 +0200
Subject: [PATCH] GEAR: add hubble flow
---
src/star_formation/GEAR/star_formation.h | 22 +++++++------------
src/star_formation/GEAR/star_formation_iact.h | 22 +++++++++++++------
2 files changed, 23 insertions(+), 21 deletions(-)
diff --git a/src/star_formation/GEAR/star_formation.h b/src/star_formation/GEAR/star_formation.h
index cfb314be43..a1d9dda2a0 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 ef42d264ee..80de642275 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 */
--
GitLab