Skip to content
GitLab
Commit e5d6d217 authored by Loic Hausammann's avatar Loic Hausammann
Browse files

GEAR: star formation fully implemented

parent bae59990
Hide whitespace changes
Inline Side-by-side
src/cooling/grackle/cooling.h
View file @ e5d6d217
......@@ -697,6 +697,17 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
xp->cooling_data.radiated_energy -= hydro_get_mass(p) * cooling_du_dt * dt;
}
/**
* @brief Compute the temperature of a #part based on the cooling function.
*
* @param phys_const #phys_const data structure.
* @param hydro_props The properties of the hydro scheme.
* @param us The internal system of units.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param p #part data.
* @param xp Pointer to the #xpart data.
*/
static INLINE float cooling_get_temperature(
const struct phys_const* restrict phys_const,
const struct hydro_props* restrict hydro_props,
......@@ -704,9 +715,30 @@ static INLINE float cooling_get_temperature(
const struct cosmology* restrict cosmo,
const struct cooling_function_data* restrict cooling,
const struct part* restrict p, const struct xpart* restrict xp) {
// TODO use the grackle library
/* Physical constants */
const double m_H = phys_const->const_proton_mass;
const double k_B = phys_const->const_boltzmann_k;
error("This function needs implementing!!");
return 0.;
/* Gas properties */
const double T_transition = hydro_props->hydrogen_ionization_temperature;
const double mu_neutral = hydro_props->mu_neutral;
const double mu_ionised = hydro_props->mu_ionised;
/* Particle temperature */
const double u = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Temperature over mean molecular weight */
const double T_over_mu = hydro_gamma_minus_one * u * m_H / k_B;
/* Are we above or below the HII -> HI transition? */
if (T_over_mu > (T_transition + 1.) / mu_ionised)
return T_over_mu * mu_ionised;
else if (T_over_mu < (T_transition - 1.) / mu_neutral)
return T_over_mu * mu_neutral;
else
return T_transition;
}
/**
......
src/star_formation/GEAR/star_formation.h
View file @ e5d6d217
......@@ -34,52 +34,12 @@
#include "star_formation_struct.h"
#include "units.h"
/**
* @brief Compute the temperature of a #part based on the cooling function.
*
* @param phys_const #phys_const data structure.
* @param hydro_props The properties of the hydro scheme.
* @param us The internal system of units.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param p #part data.
* @param xp Pointer to the #xpart data.
*/
INLINE static float get_temperature(
const struct phys_const* restrict phys_const,
const struct hydro_props* restrict hydro_props,
const struct unit_system* restrict us,
const struct cosmology* restrict cosmo,
const struct cooling_function_data* restrict cooling,
const struct part* restrict p, const struct xpart* restrict xp) {
/* Physical constants */
const double m_H = phys_const->const_proton_mass;
const double k_B = phys_const->const_boltzmann_k;
/* Gas properties */
const double T_transition = hydro_props->hydrogen_ionization_temperature;
const double mu_neutral = hydro_props->mu_neutral;
const double mu_ionised = hydro_props->mu_ionised;
/* Particle temperature */
const double u = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Temperature over mean molecular weight */
const double T_over_mu = hydro_gamma_minus_one * u * m_H / k_B;
/* Are we above or below the HII -> HI transition? */
if (T_over_mu > (T_transition + 1.) / mu_ionised)
return T_over_mu * mu_ionised;
else if (T_over_mu < (T_transition - 1.) / mu_neutral)
return T_over_mu * mu_neutral;
else
return T_transition;
}
/**
* @brief Calculate if the gas has the potential of becoming
* a star.
*
* Use the star formation criterion given by eq. 3 in Revaz & Jablonka 2018.
*
* @param starform the star formation law properties to use.
* @param p the gas particles.
* @param xp the additional properties of the gas particles.
......@@ -98,50 +58,47 @@ INLINE static int star_formation_is_star_forming(
const struct unit_system* restrict us,
const struct cooling_function_data* restrict cooling,
const struct entropy_floor_properties* restrict entropy_floor) {
/* Some useful constants */
const double G = phys_const->const_newton_G;
const double kb = phys_const->const_boltzmann_k;
const double mH = phys_const->const_proton_mass;
/* We first check if we are supposed to include turbulence estimation
* otherewise we keep 0 */
float sigma2 = 0.f;
if (starform->with_sigma > 0) {
sigma2 = xp->sf_data.sigma2;
}
/* Compute the temperature */
double const T =
get_temperature(phys_const, hydro_props, us, cosmo, cooling, p, xp);
/* Other useful values */
const int N = starform->Njeans;
const double h = p->h;
/* We suppose that the gas is neutral */
const double mu = hydro_props->mu_neutral;
/* Maximum temperature allowed (temperature criterion) */
const double T0 = starform->Max_temperature;
/* Compute density */
const double physical_density = hydro_get_physical_density(p, cosmo);
/* We compute the minimal density for star formation (see Revaz & Jablonka,
* 2018 eq (3)) */
const double rho_sfr = M_PI * (hydro_gamma * kb * T / mu / mH + sigma2) / h /
h / 4. / G / pow(N, 2. / 3.);
/* Temperature criterion for star formation eligibility */
const float temperature =
cooling_get_temperature(phys_const, hydro_props, us, cosmo, cooling, p, xp);
const float temperature_max = starform->Max_temperature;
/* Check the temperature criterion */
if (T > T0) {
return 0;
}
/* Density criterion */
else {
if (physical_density > rho_sfr) {
/* The particle is eligible (satisfy both conditions) */
return 1;
} else {
return 0;
}
/* Get the required variables */
const float G = phys_const->const_newton_G;
const float kb = phys_const->const_boltzmann_k;
const float mH = phys_const->const_proton_mass;
const float sigma2 = xp->sf_data.sigma2;
const int n_jeans_2_3 = starform->n_jeans_2_3;
const float h = p->h;
const float density = hydro_get_physical_density(p, cosmo);
// TODO use GRACKLE */
const float mu = hydro_props->mu_neutral;
/* Compute the density criterion */
const float coef = M_PI_4 / (G * n_jeans_2_3 * h * h);
const float density_criterion = coef * (hydro_gamma * kb * T / (mu * mH) + sigma2);
/* Check the density criterion */
if (density > density_criterion) {
return 1;
} else {
return 0;
}
}
/**
* @brief Compute the star-formation rate of a given particle and store
* it into the #xpart.
* @brief Compute the star-formation rate of a given particle.
*
* Nothing to do here. Everything is done in #star_formation_should_convert_to_star.
*
* @param p #part.
* @param xp the #xpart.
......@@ -153,29 +110,13 @@ INLINE static int star_formation_is_star_forming(
INLINE static void star_formation_compute_SFR(
struct part* restrict p, struct xpart* restrict xp,
const struct star_formation* starform, const struct phys_const* phys_const,
const struct cosmology* cosmo, const double dt_star) {
if (dt_star == 0.) {
xp->sf_data.SFR = 0.;
} else {
const double G = phys_const->const_newton_G;
const double c_star = starform->star_formation_rate;
const double physical_density = hydro_get_physical_density(p, cosmo);
if (physical_density != 0) {
/* We compute the star formation rate and store it (see Revaz & Jablonka,
* 2012, eq. (5)) */
/* Units are mass/time/distance³ */
xp->sf_data.SFR = c_star * physical_density *
sqrt(physical_density * 32.f * G) / sqrt(3 * M_PI);
} else {
xp->sf_data.SFR = 0.;
}
}
return;
}
const struct cosmology* cosmo, const double dt_star) {}
/**
* @brief Decides whether a particle should be converted into a
* star or not.
* Compute the star formation rate from eq. 4 in Revaz & Jablonka 2012.
*
* @param p The #part.
* @param xp The #xpart.
......@@ -187,16 +128,30 @@ INLINE static void star_formation_compute_SFR(
INLINE static int star_formation_should_convert_to_star(
struct part* p, struct xpart* xp, const struct star_formation* starform,
const struct engine* e, const double dt_star) {
/* Calculate the propability of forming a star, see Revaz & Jablonka (2012),
* eq. (4) */
const double prob = 1. - exp(xp->sf_data.SFR * dt_star * (-1.) /
hydro_get_physical_density(p, e->cosmology));
/* Get a unique random number between 0 and 1 for star formation */
const double random_number =
random_unit_interval(p->id, e->ti_current, random_number_star_formation);
/* Check that we are running a full time step */
if (dt_star == 0.) {
return 0;
}
/* Get a few variables */
const float G = phys_const->const_newton_G;
const float c_star = starform->star_formation_rate;
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_efficency * 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;
}
}
......@@ -212,19 +167,7 @@ INLINE static int star_formation_should_convert_to_star(
*/
INLINE static void star_formation_update_part_not_SFR(
struct part* p, struct xpart* xp, const struct engine* e,
const struct star_formation* starform, const int with_cosmology) {
/* Check if it is the first time steps after star formation */
if (xp->sf_data.SFR > 0.) {
/* Record the current time as an indicator of when this particle was last
star-forming. */
if (with_cosmology) {
xp->sf_data.SFR = -(e->cosmology->a);
} else {
xp->sf_data.SFR = -(e->time);
}
}
}
const struct star_formation* starform, const int with_cosmology) {}
/**
* @brief Copies the properties of the gas particle over to the
......@@ -257,6 +200,8 @@ INLINE static void star_formation_copy_properties(
} else {
sp->birth_time = e->time;
}
// TODO copy only metals
/* Store the chemistry struct in the star particle */
sp->chemistry_data = p->chemistry_data;
......@@ -265,6 +210,7 @@ INLINE static void star_formation_copy_properties(
/* Store the birth density in the star particle */
sp->birth_density = hydro_get_physical_density(p, cosmo);
/* Store the birth temperature*/
sp->birth_temperature =
get_temperature(starform->phys_const, starform->hydro_props, starform->us,
......@@ -284,24 +230,24 @@ INLINE static void starformation_init_backend(
struct swift_params* parameter_file, const struct phys_const* phys_const,
const struct unit_system* us, const struct hydro_props* hydro_props,
struct star_formation* starform) {
/* Jeans Number used in starformation elligibility criterion */
starform->Njeans =
parser_get_param_int(parameter_file, "GEARStarFormation:jeans_number");
/* Starformation efficiency (used in SFR calculation) */
starform->star_formation_rate = parser_get_param_double(
// TODO move into pressure floor
starform->n_jeans_2_3 =
parser_get_param_float(parameter_file, "GEARStarFormation:NJeans");
starform->n_jeans_2_3 = pow(starform->n_jeans_2_3, 2./3.);
/* Star formation efficiency */
starform->star_formation_efficiency = parser_get_param_double(
parameter_file, "GEARStarFormation:star_formation_efficiency");
/* Wheter we take into account local turbulence */
starform->with_sigma = parser_get_param_int(
parameter_file, "GEARStarFormation:with_turbulence_estimation");
/* Maximum temperature for starformation */
starform->Max_temperature =
/* Maximum temperature for star formation */
starform->maximal_temperature =
parser_get_param_double(parameter_file,
"GEARStarFormation:Max_temperature") *
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
;
starform->hydro_props = hydro_props;
starform->us = us;
starform->phys_const = phys_const;
"GEARStarFormation:maximal_temperature");
/* Apply unit change */
starform->maximal_temperature *=
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
}
/**
......@@ -325,9 +271,12 @@ INLINE static void starformation_print_backend(
__attribute__((always_inline)) INLINE static void star_formation_end_density(
struct part* restrict p, struct xpart* restrict xp,
const struct star_formation* cd, const struct cosmology* cosmo) {
// TODO move into pressure floor
/* To finish the turbulence estimation we devide by the density */
xp->sf_data.sigma2 /= hydro_get_physical_density(p, cosmo);
xp->sf_data.sigma2 /= pow_dimension(p->h) * hydro_get_physical_density(p, cosmo);
}
/**
* @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
*
......@@ -341,6 +290,8 @@ star_formation_part_has_no_neighbours(struct part* restrict p,
struct xpart* restrict xp,
const struct star_formation* cd,
const struct cosmology* cosmo) {
// TODO move into pressure floor
/* If part has 0 neighbours, the estimation of turbulence is 0 */
xp->sf_data.sigma2 = 0.f;
}
......@@ -361,7 +312,8 @@ star_formation_first_init_part(const struct phys_const* restrict phys_const,
const struct star_formation* data,
struct part* restrict p) {
xp->sf_data.sigma2 = 0.f;
/* Nothing special here */
star_formation_init_part(p, xp, data);
}
/**
......
src/star_formation/GEAR/star_formation_iact.h
View file @ e5d6d217
......@@ -29,6 +29,8 @@
* @brief do star_formation computation after the runner_iact_density (symmetric
* version)
*
* Compute the velocity dispersion follow eq. 2 in Revaz & Jablonka 2018.
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param hi Comoving smoothing-length of particle i.
......@@ -42,21 +44,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_star_formation(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj, struct xpart *restrict xpi,
struct xpart *restrict xpj, float a, float H) {
/* The goal here is to estimate the local turbulence */
/* Value used to evaluate the SPH kernel */
float wi;
float wj;
/* Evaluation of the SPH kernel */
kernel_eval(sqrt(r2) / hi, &wi);
kernel_eval(sqrt(r2) / hj, &wj);
/* Square of the velocity norm between particles i and j */
float norm_v2 = pow(pi->v[0] - pj->v[0], 2) + pow(pi->v[1] - pj->v[1], 2) +
pow(pi->v[2] - pj->v[2], 2);
/* Estimation of local turbulence for pi and pj, see Revaz & Jablonka, eq (2)
*/
xpi->sf_data.sigma2 += pow(hi, -3) * norm_v2 * wi * hydro_get_mass(pj);
xpj->sf_data.sigma2 += pow(hj, -3) * norm_v2 * wj * hydro_get_mass(pi);
/* Delta v */
float dv[3] = {
pi->v[0] - pj->v[0],
pi->v[1] - pj->v[1],
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];
/* Compute the velocity dispersion */
xpi->sf_data.sigma2 += norm_v2 * wi * hydro_get_mass(pj);
xpj->sf_data.sigma2 += norm_v2 * wj * hydro_get_mass(pi);
}
/**
......@@ -78,17 +85,23 @@ runner_iact_nonsym_star_formation(float r2, const float *dx, float hi, float hj,
const struct part *restrict pj,
struct xpart *restrict xpi,
const struct xpart *restrict xpj, float a,
float H)
/* The goal here is to estimate the local turbulence */
/* Value used to evaluate the SPH kernel */
{
float H) {
float wi;
/* Evaluation of the SPH kernel */
kernel_eval(sqrt(r2) / hi, &wi);
/* Square of the velocity norm */
float norm_v2 = pow(pi->v[0] - pj->v[0], 2) + pow(pi->v[1] - pj->v[1], 2) +
pow(pi->v[2] - pj->v[2], 2);
/* Estimation of local turbulence for pi */
xpi->sf_data.sigma2 += pow(hi, -3) * norm_v2 * wi * hydro_get_mass(pj);
/* Delta v */
float dv[3] = {
pi->v[0] - pj->v[0],
pi->v[1] - pj->v[1],
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];
/* Compute the velocity dispersion */
xpi->sf_data.sigma2 += norm_v2 * wi * hydro_get_mass(pj);
}
#endif /* SWIFT_GEAR_STAR_FORMATION_IACT_H */
src/star_formation/GEAR/star_formation_logger.h
View file @ e5d6d217
......@@ -195,11 +195,11 @@ INLINE static void star_formation_logger_init_log_file(
/**
* @brief Add the SFR tracer to the total active SFR of this cell
*
* Nothing to do here
*
* @param p the #part
* @param xp the #xpart
*
* Nothing to do here
*
* @param sf the SFH logger struct
* @param dt_star The length of the time-step in physical internal units.
*/
......
src/star_formation/GEAR/star_formation_struct.h
View file @ e5d6d217
......@@ -24,31 +24,26 @@
* data.
*/
struct star_formation_xpart_data {
/*!star formation rate (mass/(time*volume))*/
double SFR;
// TODO move it to the pressure floor
/*! Estimation of local turbulence (squared) */
float sigma2;
};
/* Starformation struct */
/**
* @brief Global star formation properties
*/
struct star_formation {
/*! Njeans number. We request that the SPH mass resolution is Njeans times
* smaller than the Jeans mass*/
int Njeans;
/*!star formation efficency. If the particle can create a star, it happens
* with this probablity*/
double star_formation_rate;
/*!do we include local turbulence*/
int with_sigma;
/*!Maximum temperature to allow star formation* (should be about 10'000 or
* 30'000 K*/
double Max_temperature;
/*!some other useful elements*/
const struct hydro_props* restrict hydro_props;
/*!units*/
const struct unit_system* restrict us;
/*! physical constants*/
const struct phys_const* phys_const;
// TODO move it to pressure floor
/*! Number of particle required to resolved the Jeans criterion (at power 2/3) */
float n_jeans_2_3;
/*! Maximal temperature for forming a star */
float maximal_temperature;
/*! Star formation efficiency */
float star_formation_efficiency;
};
#endif /* SWIFT_GEAR_STAR_FORMATION_STRUCT_H */
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment