diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index 996e35e41a3a6090a11ecf30a9decb391865ed63..6aeeaad3e41498b48e215f0a761ecdffb486a118 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -338,6 +338,11 @@ EAGLEChemistry: # Parameters related to star formation models ----------------------------------------------- +# GEAR star formation model (Revaz and Jablonka 2018) +GEARStarFormation: + star_formation_efficiency: 0.01 # star formation efficiency (c_*) + maximal_temperature: 3e4 # Upper limit to the temperature of a star forming particle + # EAGLE star formation model (Schaye and Dalla Vecchia 2008) EAGLEStarFormation: EOS_density_norm_H_p_cm3: 0.1 # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3. diff --git a/src/cell.c b/src/cell.c index f4903f024b10c7cc7867291b597a3bdf42997c27..0649d31b38279f3a75ec2b439f71de105425ae27 100644 --- a/src/cell.c +++ b/src/cell.c @@ -4413,7 +4413,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) { if (part_is_active(p, e)) { hydro_init_part(p, &e->s->hs); chemistry_init_part(p, e->chemistry); - star_formation_init_part(p, e->star_formation); + star_formation_init_part(p, xp, e->star_formation); tracers_after_init(p, xp, e->internal_units, e->physical_constants, with_cosmology, e->cosmology, e->hydro_properties, e->cooling_func, e->time); diff --git a/src/collectgroup.c b/src/collectgroup.c index d17608d67a291ae412728b3fcd9ea80c192802bf..8dd9e974f41c979f4c50e84fb555f9b80db19e25 100644 --- a/src/collectgroup.c +++ b/src/collectgroup.c @@ -122,7 +122,8 @@ void collectgroup1_apply(const struct collectgroup1 *grp1, struct engine *e) { e->total_nr_cells = grp1->total_nr_cells; e->total_nr_tasks = grp1->total_nr_tasks; e->tasks_per_cell_max = grp1->tasks_per_cell_max; - e->sfh = grp1->sfh; + + star_formation_logger_add_to_accumulator(&e->sfh, &grp1->sfh); } /** diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h index 1abbe5d3827726bde34502cd2b6a1d18cd309950..57f4b146c451557365e08d260ddd45276844af9c 100644 --- a/src/cooling/grackle/cooling.h +++ b/src/cooling/grackle/cooling.h @@ -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; } /** diff --git a/src/engine.c b/src/engine.c index 8f14cfd0088cff9d47bee38662a5aeff357dd8e0..33f94e9c16475b5281cdb5948e6b31e1026c7d46 100644 --- a/src/engine.c +++ b/src/engine.c @@ -5062,6 +5062,11 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params, parser_get_opt_param_double(params, "FOF:delta_time", -1.); } + /* Initialize the star formation history structure */ + if (e->policy & engine_policy_star_formation) { + star_formation_logger_accumulator_init(&e->sfh); + } + engine_init_output_lists(e, params); } diff --git a/src/engine.h b/src/engine.h index 8296933561b5f7664904d439ec6ca20f1e37a4b8..967a430871648ac4cb91f390f8b85675a314d3a8 100644 --- a/src/engine.h +++ b/src/engine.h @@ -237,7 +237,7 @@ struct engine { long long b_updates_since_rebuild; /* Star formation logger information */ - struct star_formation_history sfh; + struct star_formation_history_accumulator sfh; /* Properties of the previous step */ int step_props; diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h index 7a8844d560561dae80c676a0b5bb72b34416d080..853d2adf17bc069434562fa96ddb881f760f6830 100644 --- a/src/hydro/Gadget2/hydro_part.h +++ b/src/hydro/Gadget2/hydro_part.h @@ -155,6 +155,9 @@ struct part { /*! Black holes information (e.g. swallowing ID) */ struct black_holes_part_data black_holes_data; + /* Additional data used by the star formation */ + struct star_formation_part_data sf_data; + /* Time-step length */ timebin_t time_bin; diff --git a/src/runner.c b/src/runner.c index 7320fd8744948632cd6a08b2434768b47d837360..644d7f6171887a3a44395767353c1de84aa513d9 100644 --- a/src/runner.c +++ b/src/runner.c @@ -2272,7 +2272,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { /* Re-initialise everything */ hydro_init_part(p, hs); chemistry_init_part(p, chemistry); - star_formation_init_part(p, star_formation); + star_formation_init_part(p, xp, star_formation); tracers_after_init(p, xp, e->internal_units, e->physical_constants, with_cosmology, e->cosmology, e->hydro_properties, e->cooling_func, e->time); diff --git a/src/runner_doiact.h b/src/runner_doiact.h index 6caa287cf726f85778ef5abdc184acaf759e8b0e..1a39c8d49f2c4234c04982e255705f06ec1c5d38 100644 --- a/src/runner_doiact.h +++ b/src/runner_doiact.h @@ -884,7 +884,6 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci, const int count_i = ci->hydro.count; struct part *restrict parts_j = ci->hydro.parts; - /* Loop over the parts in ci. */ for (int pid = 0; pid < count; pid++) { diff --git a/src/space.c b/src/space.c index 49c714743728a3b93f631434d87de69a61dd716f..250af5efafa506408c9fd91c1adffb7cd96a1a21 100644 --- a/src/space.c +++ b/src/space.c @@ -4365,7 +4365,7 @@ void space_init_parts_mapper(void *restrict map_data, int count, for (int k = 0; k < count; k++) { hydro_init_part(&parts[k], hs); chemistry_init_part(&parts[k], e->chemistry); - star_formation_init_part(&parts[k], e->star_formation); + star_formation_init_part(&parts[k], &xparts[k], e->star_formation); tracers_after_init(&parts[k], &xparts[k], e->internal_units, e->physical_constants, with_cosmology, e->cosmology, e->hydro_properties, e->cooling_func, e->time); diff --git a/src/star_formation.c b/src/star_formation.c index 698a64cc636dd79f00feac3f6cc88bf519fe09c1..60cff1e2e68feaf7e71705b5079294ec478fad42 100644 --- a/src/star_formation.c +++ b/src/star_formation.c @@ -24,6 +24,7 @@ #include "part.h" #include "restart.h" #include "star_formation.h" +#include "star_formation_io.h" #include "units.h" /** diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h index 694cac1512d91ef3d640185cbf9b614773263a75..f2c77e036842b4ca040c58a6bcad1513b03a42bd 100644 --- a/src/star_formation/EAGLE/star_formation.h +++ b/src/star_formation/EAGLE/star_formation.h @@ -691,6 +691,7 @@ star_formation_first_init_part(const struct phys_const* restrict phys_const, * @param data The global star_formation information. */ __attribute__((always_inline)) INLINE static void star_formation_init_part( - struct part* restrict p, const struct star_formation* data) {} + struct part* restrict p, struct xpart* restrict xp, + const struct star_formation* data) {} #endif /* SWIFT_EAGLE_STAR_FORMATION_H */ diff --git a/src/star_formation/EAGLE/star_formation_logger.h b/src/star_formation/EAGLE/star_formation_logger.h index d634c876e52d45588ed0b93c0afc09731317c037..a843f63ce518aa9be6b1d389d42cf71c8608b903 100644 --- a/src/star_formation/EAGLE/star_formation_logger.h +++ b/src/star_formation/EAGLE/star_formation_logger.h @@ -87,6 +87,28 @@ INLINE static void star_formation_logger_add( sf_update->SFR_inactive += sf_add->SFR_inactive; } +/** + * @brief add a star formation history struct to the engine star formation + * history accumulator struct + * + * @param sf_add the star formation accumulator struct which we want to add to + * the star formation history + * @param sf_update the star formation structure which we want to update + */ +INLINE static void star_formation_logger_add_to_accumulator( + struct star_formation_history_accumulator *sf_update, + const struct star_formation_history *sf_add) { + + /* Update the SFH structure */ + sf_update->new_stellar_mass += sf_add->new_stellar_mass; + + sf_update->SFR_active += sf_add->SFR_active; + + sf_update->SFRdt_active += sf_add->SFRdt_active; + + sf_update->SFR_inactive += sf_add->SFR_inactive; +} + /** * @brief Initialize the star formation history structure in the #engine * @@ -105,6 +127,24 @@ INLINE static void star_formation_logger_init( sfh->SFR_inactive = 0.f; } +/** + * @brief Initialize the star formation history structure in the #engine + * + * @param sfh The pointer to the star formation history structure + */ +INLINE static void star_formation_logger_accumulator_init( + struct star_formation_history_accumulator *sfh) { + + /* Initialize the collecting SFH structure to zero */ + sfh->new_stellar_mass = 0.f; + + sfh->SFR_active = 0.f; + + sfh->SFRdt_active = 0.f; + + sfh->SFR_inactive = 0.f; +} + /** * @brief Write the final SFH to a file * @@ -117,7 +157,7 @@ INLINE static void star_formation_logger_init( */ INLINE static void star_formation_logger_write_to_log_file( FILE *fp, const double time, const double a, const double z, - const struct star_formation_history sf, const int step) { + const struct star_formation_history_accumulator sf, const int step) { /* Calculate the total SFR */ const float totalSFR = sf.SFR_active + sf.SFR_inactive; diff --git a/src/star_formation/EAGLE/star_formation_logger_struct.h b/src/star_formation/EAGLE/star_formation_logger_struct.h index 2a23659c4d931735d1b82a6143b3d9f871f7137a..c03a00c97ead46f552350a43574c5bbe7ac6df1b 100644 --- a/src/star_formation/EAGLE/star_formation_logger_struct.h +++ b/src/star_formation/EAGLE/star_formation_logger_struct.h @@ -34,4 +34,21 @@ struct star_formation_history { float SFRdt_active; }; +/* Starformation history struct for the engine. + Allows to integrate in time some values. + Nothing to do in EAGLE => copy of star_formation_history */ +struct star_formation_history_accumulator { + /*! Total new stellar mass */ + float new_stellar_mass; + + /*! SFR of all particles */ + float SFR_inactive; + + /*! SFR of active particles */ + float SFR_active; + + /*! SFR*dt of active particles */ + float SFRdt_active; +}; + #endif /* SWIFT_EAGLE_STAR_FORMATION_LOGGER_STRUCT_H */ diff --git a/src/star_formation/EAGLE/star_formation_struct.h b/src/star_formation/EAGLE/star_formation_struct.h index 41247e160a3eddbc9184c59b67cfa2a1d7259a05..8caac49d4b57652c5db9ae93e3789dc690e6d23f 100644 --- a/src/star_formation/EAGLE/star_formation_struct.h +++ b/src/star_formation/EAGLE/star_formation_struct.h @@ -29,4 +29,6 @@ struct star_formation_xpart_data { float SFR; }; +struct star_formation_part_data {}; + #endif /* SWIFT_EAGLE_STAR_FORMATION_STRUCT_H */ diff --git a/src/star_formation/GEAR/star_formation.h b/src/star_formation/GEAR/star_formation.h index c479feb5c66328f9fab8bf62593ca66b6658b79e..ac423a51865609460e870f65b1eeeb266182e2ef 100644 --- a/src/star_formation/GEAR/star_formation.h +++ b/src/star_formation/GEAR/star_formation.h @@ -1,6 +1,7 @@ /******************************************************************************* * This file is part of SWIFT. - * Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl) + * Coypright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) + * 2019 Fabien Jeanquartier (fabien.jeanquartier@epfl.ch) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published @@ -20,20 +21,24 @@ #define SWIFT_GEAR_STAR_FORMATION_H /* Local includes */ +#include "cooling.h" #include "cosmology.h" +#include "engine.h" #include "entropy_floor.h" #include "error.h" #include "hydro_properties.h" #include "parser.h" #include "part.h" #include "physical_constants.h" +#include "random.h" +#include "star_formation_struct.h" #include "units.h" /** * @brief Calculate if the gas has the potential of becoming * a star. * - * No star formation should occur, so return 0. + * 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. @@ -46,7 +51,7 @@ * */ INLINE static int star_formation_is_star_forming( - const struct part* restrict p, const struct xpart* restrict xp, + struct part* restrict p, struct xpart* restrict xp, const struct star_formation* starform, const struct phys_const* phys_const, const struct cosmology* cosmo, const struct hydro_props* restrict hydro_props, @@ -54,14 +59,43 @@ INLINE static int star_formation_is_star_forming( const struct cooling_function_data* restrict cooling, const struct entropy_floor_properties* restrict entropy_floor) { - return 0; + const float temperature = cooling_get_temperature(phys_const, hydro_props, us, + cosmo, cooling, p, xp); + + const float temperature_max = starform->maximal_temperature; + + /* Check the temperature criterion */ + if (temperature > temperature_max) { + return 0; + } + + /* Get the required variables */ + const float sigma2 = p->sf_data.sigma2; + const float 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 / (phys_const->const_newton_G * n_jeans_2_3 * h * h); + const float density_criterion = + coef * (hydro_gamma * phys_const->const_boltzmann_k * temperature / + (mu * phys_const->const_proton_mass) + + sigma2); + + /* Check the density criterion */ + return density > density_criterion; } /** - * @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. + * Nothing to do here. Everything is done in + * #star_formation_should_convert_to_star. * * @param p #part. * @param xp the #xpart. @@ -71,15 +105,15 @@ INLINE static int star_formation_is_star_forming( * @param dt_star The time-step of this particle. */ INLINE static void star_formation_compute_SFR( - const struct part* restrict p, struct xpart* restrict xp, + 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) {} /** * @brief Decides whether a particle should be converted into a * star or not. - * - * No SF should occur, so return 0. + + * Compute the star formation rate from eq. 4 in Revaz & Jablonka 2012. * * @param p The #part. * @param xp The #xpart. @@ -89,18 +123,38 @@ INLINE static void star_formation_compute_SFR( * @return 1 if a conversion should be done, 0 otherwise. */ INLINE static int star_formation_should_convert_to_star( - const struct part* p, const struct xpart* xp, - const struct star_formation* starform, const struct engine* e, - const double dt_star) { + struct part* p, struct xpart* xp, const struct star_formation* starform, + const struct engine* e, const double dt_star) { + + const struct phys_const* phys_const = e->physical_constants; + const struct cosmology* cosmo = e->cosmology; + + /* 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 density = hydro_get_physical_density(p, cosmo); + + /* Compute the probability */ + 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); - return 0; + /* Can we form a star? */ + return random_number < prob; } /** * @brief Update the SF properties of a particle that is not star forming. * - * Nothing to do here. - * * @param p The #part. * @param xp The #xpart. * @param e The #engine. @@ -115,8 +169,6 @@ INLINE static void star_formation_update_part_not_SFR( * @brief Copies the properties of the gas particle over to the * star particle. * - * Nothing to do here. - * * @param e The #engine * @param p the gas particles. * @param xp the additional properties of the gas particles. @@ -133,21 +185,33 @@ INLINE static void star_formation_copy_properties( const struct phys_const* phys_const, 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) { -/** - * @brief initialization of the star formation law - * - * @param parameter_file The parsed parameter file - * @param phys_const Physical constants in internal units - * @param us The current internal system of units - * @param starform the star formation law properties to initialize - * - */ -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, - const struct star_formation* starform) {} + /* Store the current mass */ + sp->mass = hydro_get_mass(p); + sp->birth.mass = sp->mass; + + /* Store either the birth_scale_factor or birth_time depending */ + if (with_cosmology) { + sp->birth_scale_factor = cosmo->a; + } else { + sp->birth_time = e->time; + } + + // TODO copy only metals + /* Store the chemistry struct in the star particle */ + sp->chemistry_data = p->chemistry_data; + + /* Store the tracers data */ + sp->tracers_data = xp->tracers_data; + + /* Store the birth density in the star particle */ + sp->birth.density = hydro_get_physical_density(p, cosmo); + + /* Store the birth temperature*/ + sp->birth.temperature = cooling_get_temperature(phys_const, hydro_props, us, + cosmo, cooling, p, xp); +} /** * @brief Prints the used parameters of the star formation law @@ -156,7 +220,6 @@ INLINE static void starformation_init_backend( */ INLINE static void starformation_print_backend( const struct star_formation* starform) { - message("Star formation law is 'GEAR'"); } @@ -164,12 +227,22 @@ INLINE static void starformation_print_backend( * @brief Finishes the density calculation. * * @param p The particle to act upon - * @param cd The global star_formation information. + * @param xp The extended particle data to act upon + * @param sf The global star_formation information. * @param cosmo The current cosmological model. */ __attribute__((always_inline)) INLINE static void star_formation_end_density( - struct part* restrict p, const struct star_formation* cd, - const struct cosmology* cosmo) {} + struct part* restrict p, const struct star_formation* sf, + const struct cosmology* cosmo) { + + // 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; +} /** * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. @@ -183,39 +256,46 @@ __attribute__((always_inline)) INLINE static void star_formation_part_has_no_neighbours(struct part* restrict p, struct xpart* restrict xp, const struct star_formation* cd, - const struct cosmology* cosmo) {} + const struct cosmology* cosmo) { + + // TODO move into pressure floor + /* If part has 0 neighbours, the estimation of turbulence is 0 */ + p->sf_data.sigma2 = 0.f; +} /** * @brief Sets the star_formation properties of the (x-)particles to a valid * start state. * - * Nothing to do here. - * + * @param p Pointer to the particle data. + * @param xp Pointer to extended particle data + * @param data The global star_formation information. + */ +__attribute__((always_inline)) INLINE static void star_formation_init_part( + struct part* restrict p, struct xpart* restrict xp, + const struct star_formation* data) { + p->sf_data.sigma2 = 0.f; +} + +/** + * @brief Sets the star_formation properties of the (x-)particles to a valid + * start state. * @param phys_const The physical constant in internal units. * @param us The unit system. * @param cosmo The current cosmological model. * @param data The global star_formation information used for this run. * @param p Pointer to the particle data. - * @param xp Pointer to the extended particle data. */ __attribute__((always_inline)) INLINE static void star_formation_first_init_part(const struct phys_const* restrict phys_const, const struct unit_system* restrict us, const struct cosmology* restrict cosmo, const struct star_formation* data, - const struct part* restrict p, - struct xpart* restrict xp) {} + struct part* restrict p, + struct xpart* restrict xp) { -/** - * @brief Sets the star_formation properties of the (x-)particles to a valid - * start state. - * - * Nothing to do here. - * - * @param p Pointer to the particle data. - * @param data The global star_formation information. - */ -__attribute__((always_inline)) INLINE static void star_formation_init_part( - struct part* restrict p, const struct star_formation* data) {} + /* Nothing special here */ + star_formation_init_part(p, xp, data); +} #endif /* SWIFT_GEAR_STAR_FORMATION_H */ diff --git a/src/star_formation/GEAR/star_formation_iact.h b/src/star_formation/GEAR/star_formation_iact.h index 749b608068650a27cbe4c9a0ca4126d2740337f3..7325b92af2840b317cf1a924a1e509b34bdffba3 100644 --- a/src/star_formation/GEAR/star_formation_iact.h +++ b/src/star_formation/GEAR/star_formation_iact.h @@ -1,6 +1,7 @@ /******************************************************************************* * This file is part of SWIFT. - * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) + * Coypright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) + * 2019 Fabien Jeanquartier (fabien.jeanquartier@epfl.ch) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published @@ -28,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. @@ -39,7 +42,28 @@ */ __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, float a, float H) {} + struct part *restrict pj, float a, float H) { + + float wi; + float wj; + /* Evaluation of the SPH kernel */ + kernel_eval(sqrt(r2) / hi, &wi); + kernel_eval(sqrt(r2) / hj, &wj); + + /* Delta v */ + float dv[3] = {pi->v[0] - pj->v[0], pi->v[1] - pj->v[1], pi->v[2] - pj->v[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 += sigma2 * wi * hydro_get_mass(pj); + pj->sf_data.sigma2 += sigma2 * wj * hydro_get_mass(pi); +} /** * @brief do star_formation computation after the runner_iact_density (non @@ -58,6 +82,23 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_star_formation(float r2, const float *dx, float hi, float hj, struct part *restrict pi, const struct part *restrict pj, float a, - float H) {} + float H) { + float wi; + /* Evaluation of the SPH kernel */ + kernel_eval(sqrt(r2) / hi, &wi); + + /* Delta v */ + float dv[3] = {pi->v[0] - pj->v[0], pi->v[1] - pj->v[1], pi->v[2] - pj->v[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 += sigma2 * wi * hydro_get_mass(pj); +} #endif /* SWIFT_GEAR_STAR_FORMATION_IACT_H */ diff --git a/src/star_formation/GEAR/star_formation_io.h b/src/star_formation/GEAR/star_formation_io.h index 6ef04c49c4abcd00175aaa164271628a9ff89360..3e021f7844c1deaeca40d7144d6f7b69cb6c2bdb 100644 --- a/src/star_formation/GEAR/star_formation_io.h +++ b/src/star_formation/GEAR/star_formation_io.h @@ -1,6 +1,7 @@ /******************************************************************************* * This file is part of SWIFT. - * Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl) + * Coypright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) + * 2019 Fabien Jeanquartier (fabien.jeanquartier@epfl.ch) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published @@ -37,8 +38,40 @@ __attribute__((always_inline)) INLINE static int star_formation_write_particles( const struct part* parts, const struct xpart* xparts, struct io_props* list) { - + /* Nothing to write here */ return 0; } +/** + * @brief initialization of the star formation law + * + * @param parameter_file The parsed parameter file + * @param phys_const Physical constants in internal units + * @param us The current internal system of units + * @param starform the star formation law properties to initialize + * + */ +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) { + + // 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"); + + /* Maximum temperature for star formation */ + starform->maximal_temperature = parser_get_param_double( + parameter_file, "GEARStarFormation:maximal_temperature"); + + /* Apply unit change */ + starform->maximal_temperature *= + units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE); +} + #endif /* SWIFT_STAR_FORMATION_GEAR_IO_H */ diff --git a/src/star_formation/GEAR/star_formation_logger.h b/src/star_formation/GEAR/star_formation_logger.h index 5b9e033d21d3d202f3c289d1dd6a843ba17fa524..84475909c0f524a4f930a48dbcc3b0943719f8b0 100644 --- a/src/star_formation/GEAR/star_formation_logger.h +++ b/src/star_formation/GEAR/star_formation_logger.h @@ -1,6 +1,7 @@ /******************************************************************************* * This file is part of SWIFT. - * Copyright (c) 2019 Folkert Nobels (nobels@strw.leidenuniv.nl) + * Coypright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) + * 2019 Fabien Jeanquartier (fabien.jeanquartier@epfl.ch) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published @@ -27,16 +28,28 @@ #include "hydro.h" #include "part.h" #include "star_formation_logger_struct.h" +#include "units.h" /** - * @brief Update the stellar mass in the current cell after creating - * the new star particle spart sp + * @brief Update the stellar quantities in the current cell after creating + * the new star particle spart sp. * + * @param time_step, the current time step of the simulation * @param sp new created star particle * @param sf the star_formation_history struct of the current cell */ INLINE static void star_formation_logger_log_new_spart( - struct spart *sp, struct star_formation_history *sf) {} + const struct spart *sp, struct star_formation_history *sf) { + + /* Add mass of created sparticle to the total stellar mass in this cell */ + sf->new_stellar_mass += sp->mass; + + /* Increase the number of stars */ + sf->number_new_stars += 1; + + /* No need to deal with the integrated quantities, only the engine's one is + * updated */ +} /** * @brief Initialize the star formation history struct in the case the cell is @@ -45,41 +58,42 @@ INLINE static void star_formation_logger_log_new_spart( * @param sf the star_formation_history struct we want to initialize */ INLINE static void star_formation_logger_log_inactive_cell( - struct star_formation_history *sf) {} + struct star_formation_history *sf) { -/** - * @brief add a star formation history struct to an other star formation history - * struct - * - * @param sf_add the star formation struct which we want to add to the star - * formation history - * @param sf_update the star formation structure which we want to update - */ -INLINE static void star_formation_logger_add( - struct star_formation_history *sf_update, - const struct star_formation_history *sf_add) {} + /* Initialize the stellar mass to zero */ + sf->new_stellar_mass = 0.f; + /* initialize number of stars to zero*/ + sf->number_new_stars = 0; +} /** * @brief Initialize the star formation history structure in the #engine * * @param sfh The pointer to the star formation history structure */ INLINE static void star_formation_logger_init( - struct star_formation_history *sfh) {} + struct star_formation_history *sfh) { + /* Initialize the collecting SFH structure to zero */ + sfh->new_stellar_mass = 0.f; + sfh->number_new_stars = 0; +} /** - * @brief Write the final SFH to a file + * @brief add a star formation history struct to an other star formation history + * struct * - * @param fp The file to write to. - * @param time the simulation time (time since Big Bang) in internal units. - * @param a the scale factor. - * @param z the redshift. - * @param sf the #star_formation_history struct. - * @param step The time-step of the simulation. + * @param sf_add the star formation struct which we want to add to the star + * formation history + * @param sf_update the star formation structure which we want to update */ -INLINE static void star_formation_logger_write_to_log_file( - FILE *fp, const double time, const double a, const double z, - const struct star_formation_history sf, const int step) {} +INLINE static void star_formation_logger_add( + struct star_formation_history *sf_update, + const struct star_formation_history *sf_add) { + + /* Update the SFH structure */ + sf_update->number_new_stars += sf_add->number_new_stars; + sf_update->new_stellar_mass += sf_add->new_stellar_mass; +} /** * @brief Initialize the SFH logger file @@ -90,13 +104,47 @@ INLINE static void star_formation_logger_write_to_log_file( */ INLINE static void star_formation_logger_init_log_file( FILE *fp, const struct unit_system *restrict us, - const struct phys_const *phys_const) {} + const struct phys_const *phys_const) { + + /* Write some general text to the logger file */ + fprintf(fp, "# Star Formation History Logger file\n"); + fprintf(fp, "######################################################\n"); + fprintf(fp, "# The quantities are all given in internal physical units!\n"); + fprintf(fp, "#\n"); + fprintf(fp, "# (0) Simulation step\n"); + fprintf(fp, + "# (1) Time since Big Bang (cosmological run), Time since start of " + "the simulation (non-cosmological run).\n"); + fprintf(fp, "# Unit = %e seconds\n", us->UnitTime_in_cgs); + fprintf(fp, "# Unit = %e yr or %e Myr\n", 1.f / phys_const->const_year, + 1.f / phys_const->const_year / 1e6); + fprintf(fp, "# (2) Scale factor (no unit)\n"); + fprintf(fp, "# (3) Redshift (no unit)\n"); + fprintf(fp, + "# (4) Total number of stars formed in the simulation (no unit)\n"); + fprintf(fp, "# (5) Total stellar mass formed in the simulation.\n"); + fprintf(fp, "# Unit = %e gram\n", us->UnitMass_in_cgs); + fprintf(fp, "# Unit = %e solar mass\n", + 1.f / phys_const->const_solar_mass); + fprintf(fp, + "# (6) Number of stars formed in the current time step (no unit).\n"); + fprintf(fp, "# (7) Mass of stars formed in the current time step.\n"); + fprintf(fp, "# Unit = %e gram\n", us->UnitMass_in_cgs); + fprintf(fp, "# Unit = %e solar mass\n", + 1.f / phys_const->const_solar_mass); + fprintf(fp, + "# (0) (1) (2) (3) (4) " + " (5) (6) (7)\n"); +} /** * @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 + * * @param sf the SFH logger struct * @param dt_star The length of the time-step in physical internal units. */ @@ -108,6 +156,8 @@ INLINE static void star_formation_logger_log_active_part( * @brief Add the SFR tracer to the total inactive SFR of this cell as long as * the SFR tracer is larger than 0 * + * Nothing to do here + * * @param p the #part * @param xp the #xpart * @param sf the SFH logger struct @@ -116,4 +166,57 @@ INLINE static void star_formation_logger_log_inactive_part( const struct part *p, const struct xpart *xp, struct star_formation_history *sf) {} +/** + * @brief add a star formation history struct to an other star formation history + * struct + * + * @param sf_add the star formation struct which we want to add to the star + * formation history + * @param sf_update the star formation structure which we want to update + */ +INLINE static void star_formation_logger_add_to_accumulator( + struct star_formation_history_accumulator *sf_update, + const struct star_formation_history *sf_add) { + + /* Update the SFH structure */ + sf_update->number_new_stars = sf_add->number_new_stars; + sf_update->new_stellar_mass = sf_add->new_stellar_mass; + sf_update->total_number_stars += sf_add->number_new_stars; + sf_update->total_stellar_mass += sf_add->new_stellar_mass; +} + +/** + * @brief Write the final SFH to a file + * + * @param fp The file to write to. + * @param time the simulation time (time since Big Bang) in internal units. + * @param a the scale factor. + * @param z the redshift. + * @param sf the #star_formation_history struct. + * @param step The time-step of the simulation. + */ +INLINE static void star_formation_logger_write_to_log_file( + FILE *fp, const double time, const double a, const double z, + struct star_formation_history_accumulator sf, const int step) { + + fprintf(fp, "%6d %16e %12.7f %14e %14ld %14e %14ld %14e\n", step, time, a, z, + sf.total_number_stars, sf.total_stellar_mass, sf.number_new_stars, + sf.new_stellar_mass); +} + +/** + * @brief Initialize the star formation history struct in the #engine when + * starting the simulation. + * + * @param sfh the star_formation_history struct we want to initialize + */ +INLINE static void star_formation_logger_accumulator_init( + struct star_formation_history_accumulator *sfh) { + /* Initialize all values to zero */ + sfh->new_stellar_mass = 0.f; + sfh->number_new_stars = 0; + sfh->total_number_stars = 0; + sfh->total_stellar_mass = 0.f; +} + #endif /* SWIFT_GEAR_STARFORMATION_LOGGER_H */ diff --git a/src/star_formation/GEAR/star_formation_logger_struct.h b/src/star_formation/GEAR/star_formation_logger_struct.h index 04b5dfdc038f7b684cfb1f2079d13eb312624b3f..63e3af06a7cd5af375662a1ba28dc5000a69dc3f 100644 --- a/src/star_formation/GEAR/star_formation_logger_struct.h +++ b/src/star_formation/GEAR/star_formation_logger_struct.h @@ -1,6 +1,7 @@ /******************************************************************************* * This file is part of SWIFT. - * Copyright (c) 2019 Folkert Nobels (nobels@strw.leidenuniv.nl) + * Coypright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) + * 2019 Fabien Jeanquartier (fabien.jeanquartier@epfl.ch) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published @@ -19,7 +20,33 @@ #ifndef SWIFT_GEAR_STAR_FORMATION_LOGGER_STRUCT_H #define SWIFT_GEAR_STAR_FORMATION_LOGGER_STRUCT_H -/* Starformation history struct */ -struct star_formation_history {}; +/** + * Structure containing the star formation information from the cells. + */ +struct star_formation_history { + /*! Stellar mass created in the current timestep */ + float new_stellar_mass; -#endif /* SWIFT_NONE_STAR_FORMATION_STRUCT_H */ + /*! Number of stars created in this timestep */ + long int number_new_stars; +}; + +/** + * Structure containing the global star formation information (including time + * integrated variables). + */ +struct star_formation_history_accumulator { + /*! Total stellar mass from the begining of the simulation */ + float total_stellar_mass; + + /*! Total number of stars */ + long int total_number_stars; + + /*! Stellar mass created in the current timestep */ + float new_stellar_mass; + + /*! Number of stars created in this timestep */ + long int number_new_stars; +}; + +#endif /* SWIFT_GEAR_STAR_FORMATION_LOGGER_STRUCT_H */ diff --git a/src/star_formation/GEAR/star_formation_struct.h b/src/star_formation/GEAR/star_formation_struct.h index 9b4e216fd2955f29d89dade6ee46c0e1af715cdb..2e9a7548f83ca6ae9bb78ee7bcf4be69a0a31489 100644 --- a/src/star_formation/GEAR/star_formation_struct.h +++ b/src/star_formation/GEAR/star_formation_struct.h @@ -25,7 +25,27 @@ */ struct star_formation_xpart_data {}; -/* Starformation struct */ -struct star_formation {}; +struct star_formation_part_data { + // TODO move it to the pressure floor + /*! Estimation of local turbulence (squared) */ + float sigma2; +}; + +/** + * @brief Global star formation properties + */ +struct star_formation { + + // 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 */ diff --git a/src/star_formation/none/star_formation.h b/src/star_formation/none/star_formation.h index 412bda2cc356729bb9f3e3657689fb5bb7b36d8f..6b83010ec5f8935778af8c9ad21010ff3452fc0e 100644 --- a/src/star_formation/none/star_formation.h +++ b/src/star_formation/none/star_formation.h @@ -220,6 +220,7 @@ star_formation_first_init_part(const struct phys_const* restrict phys_const, * @param data The global star_formation information. */ __attribute__((always_inline)) INLINE static void star_formation_init_part( - struct part* restrict p, const struct star_formation* data) {} + struct part* restrict p, struct xpart* restrict xp, + const struct star_formation* data) {} #endif /* SWIFT_NONE_STAR_FORMATION_H */ diff --git a/src/star_formation/none/star_formation_logger.h b/src/star_formation/none/star_formation_logger.h index b4e6987c03d295348fc8c22d66cb20d10e54378c..552df0c6cae533d1eb2678cb23dd675e0a058715 100644 --- a/src/star_formation/none/star_formation_logger.h +++ b/src/star_formation/none/star_formation_logger.h @@ -59,6 +59,18 @@ INLINE static void star_formation_logger_add( struct star_formation_history *sf_update, const struct star_formation_history *sf_add) {} +/** + * @brief add a star formation history accumulator struct to an other star + * formation history struct + * + * @param sf_add the star formation accumulator struct which we want to add to + * the star formation history + * @param sf_update the star formation structure which we want to update + */ +INLINE static void star_formation_logger_add_to_accumulator( + struct star_formation_history_accumulator *sf_update, + const struct star_formation_history *sf_add) {} + /** * @brief Initialize the star formation history structure in the #engine * @@ -67,6 +79,14 @@ INLINE static void star_formation_logger_add( INLINE static void star_formation_logger_init( struct star_formation_history *sfh) {} +/** + * @brief Initialize the star formation history structure in the #engine + * + * @param sfh The pointer to the star formation history structure + */ +INLINE static void star_formation_logger_accumulator_init( + struct star_formation_history_accumulator *sfh) {} + /** * @brief Write the final SFH to a file * @@ -74,12 +94,12 @@ INLINE static void star_formation_logger_init( * @param time the simulation time (time since Big Bang) in internal units. * @param a the scale factor. * @param z the redshift. - * @param sf the #star_formation_history struct. + * @param sf the #star_formation_history_accumulator struct. * @param step The time-step of the simulation. */ INLINE static void star_formation_logger_write_to_log_file( FILE *fp, const double time, const double a, const double z, - const struct star_formation_history sf, const int step) {} + const struct star_formation_history_accumulator sf, const int step) {} /** * @brief Initialize the SFH logger file diff --git a/src/star_formation/none/star_formation_logger_struct.h b/src/star_formation/none/star_formation_logger_struct.h index 9efda271da96faf2088169fd75d0e3c01247a429..b60c64f2eb47894db828a1bde0ca20803892c7db 100644 --- a/src/star_formation/none/star_formation_logger_struct.h +++ b/src/star_formation/none/star_formation_logger_struct.h @@ -22,4 +22,10 @@ /* Starformation history struct */ struct star_formation_history {}; +/* Starformation history accumulator struct. + This structure is only defined in the engine and + allows the user to integrate some quantities over time. +*/ +struct star_formation_history_accumulator {}; + #endif /* SWIFT_NONE_STAR_FORMATION_STRUCT_H */ diff --git a/src/star_formation/none/star_formation_struct.h b/src/star_formation/none/star_formation_struct.h index 27a2adaf83d0a02a0d08e7eef8b45bea630689e4..2f5241a58caf1ca70fa98a40d467c8ff5a3237f7 100644 --- a/src/star_formation/none/star_formation_struct.h +++ b/src/star_formation/none/star_formation_struct.h @@ -25,4 +25,10 @@ */ struct star_formation_xpart_data {}; +/** + * @brief Star-formation-related properties stored in the particle + * data. + */ +struct star_formation_part_data {}; + #endif /* SWIFT_NONE_STAR_FORMATION_STRUCT_H */ diff --git a/src/stars.h b/src/stars.h index dd8390e0206580fc2a07a08e51bb69c6ee5ab5ed..dea6e07a87cabd7d1778ec2a850be4f3b16b04b0 100644 --- a/src/stars.h +++ b/src/stars.h @@ -29,6 +29,9 @@ #elif defined(STARS_EAGLE) #include "./stars/EAGLE/stars.h" #include "./stars/EAGLE/stars_iact.h" +#elif defined(STARS_GEAR) +#include "./stars/GEAR/stars.h" +#include "./stars/GEAR/stars_iact.h" #else #error "Invalid choice of star model" #endif