diff --git a/src/runner.c b/src/runner.c index f650f8f459c13a3fad7902b059feab736ad9ef48..099817e0107dbcb16cb2edbbeb8474d965d032c4 100644 --- a/src/runner.c +++ b/src/runner.c @@ -475,6 +475,8 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { const struct hydro_props *restrict hydro_props = e->hydro_properties; const struct unit_system *restrict us = e->internal_units; struct cooling_function_data *restrict cooling = e->cooling_func; + const double time_base = e->time_base; + const integertime_t ti_current = e->ti_current; TIMER_TIC; @@ -496,12 +498,27 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { if (part_is_active(p, e)) { + /* Calculate the time step of the current particle*/ + double dt_star; + if (with_cosmology) { + const integertime_t ti_step = get_integer_timestep(p->time_bin); + const integertime_t ti_begin = + get_integer_time_begin(ti_current - 1, p->time_bin); + + dt_star = + cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step); + + } else { + dt_star = get_timestep(p->time_bin, time_base); + } + // const float rho = hydro_get_physical_density(p, cosmo); if (star_formation_convert_to_star(e, starform, p, xp, constants, cosmo, - hydro_props, us, cooling)) { + hydro_props, us, cooling, dt_star)) { /* Convert your particle to a star */ struct spart* sp = cell_convert_part_to_spart(e, c, p, xp); + /* Copy the properties of the gas particle to the star particle */ star_formation_copy_properties(e, p, xp, sp, starform, constants, cosmo, with_cosmology); // struct spart *sp = cell_conert_part_to_spart(c, p, ...); diff --git a/src/starformation/schaye08/starformation.h b/src/starformation/schaye08/starformation.h index 344d822799f3605b543dfcb1de06d4430ae7cab9..cb5e546de1d24bd6ef891220847440561eeb4adb 100644 --- a/src/starformation/schaye08/starformation.h +++ b/src/starformation/schaye08/starformation.h @@ -156,7 +156,6 @@ INLINE static int star_formation_potential_to_become_star( if (particle_density < rho_crit_times_min_over_den) return 0; - /* In this case there are actually multiple possibilities * because we also need to check if the physical density exceeded * the appropriate limit */ @@ -165,14 +164,14 @@ INLINE static int star_formation_potential_to_become_star( double density_threshold_metal_dep = starform->density_threshold * pow(Z * starform->Z0_inv, starform->n_Z0); - double density_threshold_current = min(density_threshold_metal_dep, starform->density_threshold_max); + /* Calculate the maximum between both and convert to mass density instead of number density*/ + double density_threshold_current = min(density_threshold_metal_dep, starform->density_threshold_max) * phys_const->const_proton_mass; /* Check if it exceeded the maximum density */ - if (particle_density < density_threshold_current) + if (particle_density*p->chemistry_data.smoothed_metal_mass_fraction[0] < density_threshold_current) return 0; - /* Calculate the temperature */ const double temperature = cooling_get_temperature(phys_const, hydro_props, us, cosmo, cooling, p, xp); @@ -200,18 +199,19 @@ INLINE static int star_formation_convert_to_star( const struct phys_const* const phys_const, const struct cosmology* cosmo, const struct hydro_props* restrict hydro_props, const struct unit_system* restrict us, - const struct cooling_function_data* restrict cooling) { + const struct cooling_function_data* restrict cooling, + const double dt_star) { if (star_formation_potential_to_become_star( starform, p, xp, phys_const, cosmo, hydro_props, us, cooling)) { /* Get the pressure */ const double pressure = starform->EOS_pressure_norm * - pow(p->rho / starform->EOS_density_norm, starform->polytropic_index); + pow(hydro_get_physical_density(p,cosmo) / starform->EOS_density_norm, starform->polytropic_index); /* Calculate the propability of forming a star */ const double prop = starform->SF_normalization * - pow(pressure, starform->SF_power_law) * p->time_bin; + pow(pressure, starform->SF_power_law) * dt_star; /* Calculate the seed */ unsigned int seed = (p->id + e->ti_current) % 8191; @@ -219,6 +219,8 @@ INLINE static int star_formation_convert_to_star( /* Generate a random number between 0 and 1. */ const double randomnumber = rand_r(&seed) * starform->inv_RAND_MAX; + message("Passed whole boundary thing! random number = %e, prop = %e time_bin %d", randomnumber, prop, p->time_bin); + /* Calculate if we form a star */ return (prop > randomnumber); @@ -304,7 +306,7 @@ INLINE static void starformation_init_backend( /* Conversion of number density from cgs */ static const float dimension_numb_den[5] = {0, -3, 0, 0, 0}; - const double conversion_numb_density = + const double conversion_numb_density = 1/ units_general_cgs_conversion_factor(us, dimension_numb_den); /* Quantities that have to do with the Normal Kennicutt-