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 1fc6531656c94c566e1ffac502b5023e09094872..55f1c8d4271b16f63ce00e6d6cac9eaee5a778a6 100644 --- a/src/star_formation/EAGLE/star_formation.h +++ b/src/star_formation/EAGLE/star_formation.h @@ -690,6 +690,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 dfe645718d689841f89cf592194d435af299a642..3275e1a4c43ce232d73dbb2f331c1c6e81118ec1 100644 --- a/src/star_formation/none/star_formation.h +++ b/src/star_formation/none/star_formation.h @@ -219,6 +219,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 diff --git a/src/stars/GEAR/stars.h b/src/stars/GEAR/stars.h new file mode 100644 index 0000000000000000000000000000000000000000..467aaa164ba9762d85f8dfae85db86ff76ae779e --- /dev/null +++ b/src/stars/GEAR/stars.h @@ -0,0 +1,172 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Coypright (c) 2019 Loic Hausammann (loic.hausammann@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 + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_GEAR_STARS_H +#define SWIFT_GEAR_STARS_H + +#include <float.h> +#include "minmax.h" + +/** + * @brief Computes the gravity time-step of a given star particle. + * + * @param sp Pointer to the s-particle data. + */ +__attribute__((always_inline)) INLINE static float stars_compute_timestep( + const struct spart* const sp) { + + return FLT_MAX; +} + +/** + * @brief Initialises the s-particles for the first time + * + * This function is called only once just after the ICs have been + * read in to do some conversions. + * + * @param sp The particle to act upon + * @param stars_properties The properties of the stellar model. + */ +__attribute__((always_inline)) INLINE static void stars_first_init_spart( + struct spart* sp, const struct stars_props* stars_properties) { + + sp->time_bin = 0; +} + +/** + * @brief Prepares a s-particle for its interactions + * + * @param sp The particle to act upon + */ +__attribute__((always_inline)) INLINE static void stars_init_spart( + struct spart* sp) { + +#ifdef DEBUG_INTERACTIONS_STARS + for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i) + sp->ids_ngbs_density[i] = -1; + sp->num_ngb_density = 0; +#endif + + sp->density.wcount = 0.f; + sp->density.wcount_dh = 0.f; +} + +/** + * @brief Predict additional particle fields forward in time when drifting + * + * @param sp The particle + * @param dt_drift The drift time-step for positions. + */ +__attribute__((always_inline)) INLINE static void stars_predict_extra( + struct spart* restrict sp, float dt_drift) {} + +/** + * @brief Sets the values to be predicted in the drifts to their values at a + * kick time + * + * @param sp The particle. + */ +__attribute__((always_inline)) INLINE static void stars_reset_predicted_values( + struct spart* restrict sp) {} + +/** + * @brief Finishes the calculation of (non-gravity) forces acting on stars + * + * @param sp The particle to act upon + */ +__attribute__((always_inline)) INLINE static void stars_end_feedback( + struct spart* sp) {} + +/** + * @brief Kick the additional variables + * + * @param sp The particle to act upon + * @param dt The time-step for this kick + */ +__attribute__((always_inline)) INLINE static void stars_kick_extra( + struct spart* sp, float dt) {} + +/** + * @brief Finishes the calculation of density on stars + * + * @param sp The particle to act upon + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void stars_end_density( + struct spart* sp, const struct cosmology* cosmo) { + + /* Some smoothing length multiples. */ + const float h = sp->h; + const float h_inv = 1.0f / h; /* 1/h */ + const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ + const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */ + + /* Finish the calculation by inserting the missing h-factors */ + sp->density.wcount *= h_inv_dim; + sp->density.wcount_dh *= h_inv_dim_plus_one; +} + +/** + * @brief Sets all particle fields to sensible values when the #spart has 0 + * ngbs. + * + * @param sp The particle to act upon + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void stars_spart_has_no_neighbours( + struct spart* restrict sp, const struct cosmology* cosmo) { + + /* Some smoothing length multiples. */ + const float h = sp->h; + const float h_inv = 1.0f / h; /* 1/h */ + const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ + + /* Re-set problematic values */ + sp->density.wcount = kernel_root * h_inv_dim; + sp->density.wcount_dh = 0.f; +} + +/** + * @brief Reset acceleration fields of a particle + * + * This is the equivalent of hydro_reset_acceleration. + * We do not compute the acceleration on star, therefore no need to use it. + * + * @param p The particle to act upon + */ +__attribute__((always_inline)) INLINE static void stars_reset_feedback( + struct spart* restrict p) { + +#ifdef DEBUG_INTERACTIONS_STARS + for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i) + p->ids_ngbs_force[i] = -1; + p->num_ngb_force = 0; +#endif +} + +/** + * @brief Initializes constants related to stellar evolution, initializes imf, + * reads and processes yield tables + * + * @param params swift_params parameters structure + * @param stars stars_props data structure + */ +inline static void stars_evolve_init(struct swift_params* params, + struct stars_props* restrict stars) {} + +#endif /* SWIFT_GEAR_STARS_H */ diff --git a/src/stars/GEAR/stars_debug.h b/src/stars/GEAR/stars_debug.h new file mode 100644 index 0000000000000000000000000000000000000000..41953367f6c8771ffd7460cee493c8330ecd874b --- /dev/null +++ b/src/stars/GEAR/stars_debug.h @@ -0,0 +1,31 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Coypright (c) 2019 Loic Hausammann (loic.hausammann@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 + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_GEAR_STARS_DEBUG_H +#define SWIFT_GEAR_STARS_DEBUG_H + +__attribute__((always_inline)) INLINE static void stars_debug_particle( + const struct spart* p) { + printf( + "x=[%.3e,%.3e,%.3e], " + "v_full=[%.3e,%.3e,%.3e] p->mass=%.3e \n t_begin=%d, t_end=%d\n", + p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2], + p->mass, p->ti_begin, p->ti_end); +} + +#endif /* SWIFT_GEAR_STARS_DEBUG_H */ diff --git a/src/stars/GEAR/stars_iact.h b/src/stars/GEAR/stars_iact.h new file mode 100644 index 0000000000000000000000000000000000000000..c7bda43fc0c66fcc890e6b05bacf871edfd10b5a --- /dev/null +++ b/src/stars/GEAR/stars_iact.h @@ -0,0 +1,94 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Coypright (c) 2019 Loic Hausammann (loic.hausammann@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 + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_GEAR_STARS_IACT_H +#define SWIFT_GEAR_STARS_IACT_H + +/** + * @brief Density interaction between two particles (non-symmetric). + * + * @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. + * @param hj Comoving smoothing-length of particle j. + * @param si First sparticle. + * @param pj Second particle (not updated). + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void +runner_iact_nonsym_stars_density(const float r2, const float *dx, + const float hi, const float hj, + struct spart *restrict si, + const struct part *restrict pj, const float a, + const float H) { + + float wi, wi_dx; + + /* Get r and 1/r. */ + const float r_inv = 1.0f / sqrtf(r2); + const float r = r2 * r_inv; + + /* Compute the kernel function */ + const float hi_inv = 1.0f / hi; + const float ui = r * hi_inv; + kernel_deval(ui, &wi, &wi_dx); + + /* Compute contribution to the number of neighbours */ + si->density.wcount += wi; + si->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx); + +#ifdef DEBUG_INTERACTIONS_STARS + /* Update ngb counters */ + if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_STARS) + si->ids_ngbs_density[si->num_ngb_density] = pj->id; + + /* Update ngb counters */ + ++si->num_ngb_density; +#endif +} + +/** + * @brief Feedback interaction between two particles (non-symmetric). + * + * @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. + * @param hj Comoving smoothing-length of particle j. + * @param si First sparticle. + * @param pj Second particle (not updated). + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void +runner_iact_nonsym_stars_feedback(const float r2, const float *dx, + const float hi, const float hj, + struct spart *restrict si, + struct part *restrict pj, const float a, + const float H) { +#ifdef DEBUG_INTERACTIONS_STARS + /* Update ngb counters */ + if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_STARS) + si->ids_ngbs_force[si->num_ngb_force] = pj->id; + + /* Update ngb counters */ + ++si->num_ngb_force; +#endif +} + +#endif diff --git a/src/stars/GEAR/stars_io.h b/src/stars/GEAR/stars_io.h new file mode 100644 index 0000000000000000000000000000000000000000..ebd72aa50a4194bf8c6f747e55d265ace0550c35 --- /dev/null +++ b/src/stars/GEAR/stars_io.h @@ -0,0 +1,248 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Coypright (c) 2019 Loic Hausammann (loic.hausammann@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 + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_GEAR_STARS_IO_H +#define SWIFT_GEAR_STARS_IO_H + +#include "io_properties.h" +#include "stars_part.h" + +/** + * @brief Specifies which s-particle fields to read from a dataset + * + * @param sparts The s-particle array. + * @param list The list of i/o properties to read. + * @param num_fields The number of i/o fields to read. + */ +INLINE static void stars_read_particles(struct spart *sparts, + struct io_props *list, + int *num_fields) { + + /* Say how much we want to read */ + *num_fields = 5; + + /* List what we want to read */ + list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY, + UNIT_CONV_LENGTH, sparts, x); + list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY, + UNIT_CONV_SPEED, sparts, v); + list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS, + sparts, mass); + list[3] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY, + UNIT_CONV_NO_UNITS, sparts, id); + list[4] = io_make_input_field("SmoothingLength", FLOAT, 1, OPTIONAL, + UNIT_CONV_LENGTH, sparts, h); +} + +/** + * @brief Specifies which s-particle fields to write to a dataset + * + * @param sparts The s-particle array. + * @param list The list of i/o properties to write. + * @param num_fields The number of i/o fields to write. + * @param with_cosmology Is it a cosmological run? + */ +INLINE static void stars_write_particles(const struct spart *sparts, + struct io_props *list, int *num_fields, + const int with_cosmology) { + + /* Say how much we want to write */ + *num_fields = 9; + + /* List what we want to write */ + list[0] = + io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1., + sparts, x, "Co-moving positions of the particles"); + + list[1] = io_make_output_field( + "Velocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, sparts, v, + "Peculiar velocities of the stars. This is (a * dx/dt) where x is the " + "co-moving positions of the particles"); + + list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f, + sparts, mass, "Masses of the particles"); + + list[3] = + io_make_output_field("ParticleIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, + sparts, id, "Unique IDs of the particles"); + + list[4] = io_make_output_field( + "SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, sparts, h, + "Co-moving smoothing lengths (FWHM of the kernel) of the particles"); + + if (with_cosmology) { + list[5] = io_make_output_field( + "BirthScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, sparts, + birth_scale_factor, "Scale-factors at which the stars were born"); + } else { + list[5] = io_make_output_field("BirthTimes", FLOAT, 1, UNIT_CONV_TIME, 0.f, + sparts, birth_time, + "Times at which the stars were born"); + } + + list[6] = io_make_output_field( + "BirthDensities", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, sparts, birth.density, + "Physical densities at the time of birth of the gas particles that " + "turned into stars (note that " + "we store the physical density at the birth redshift, no conversion is " + "needed)"); + + list[7] = + io_make_output_field("BirthTemperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, + 0.f, sparts, birth.temperature, + "Temperatures at the time of birth of the gas " + "particles that turned into stars"); + + list[7] = io_make_output_field("BirthMasses", FLOAT, 1, UNIT_CONV_MASS, 0.f, + sparts, birth.mass, + "Masses of the star particles at birth time"); + +#ifdef DEBUG_INTERACTIONS_STARS + + list += *num_fields; + *num_fields += 4; + + list[0] = io_make_output_field("Num_ngb_density", INT, 1, UNIT_CONV_NO_UNITS, + sparts, num_ngb_density); + list[1] = io_make_output_field("Num_ngb_force", INT, 1, UNIT_CONV_NO_UNITS, + sparts, num_ngb_force); + list[2] = io_make_output_field("Ids_ngb_density", LONGLONG, + MAX_NUM_OF_NEIGHBOURS_STARS, + UNIT_CONV_NO_UNITS, sparts, ids_ngbs_density); + list[3] = io_make_output_field("Ids_ngb_force", LONGLONG, + MAX_NUM_OF_NEIGHBOURS_STARS, + UNIT_CONV_NO_UNITS, sparts, ids_ngbs_force); +#endif +} + +/** + * @brief Initialize the global properties of the stars scheme. + * + * By default, takes the values provided by the hydro. + * + * @param sp The #stars_props. + * @param phys_const The physical constants in the internal unit system. + * @param us The internal unit system. + * @param params The parsed parameters. + * @param p The already read-in properties of the hydro scheme. + * @param cosmo The cosmological model. + */ +INLINE static void stars_props_init(struct stars_props *sp, + const struct phys_const *phys_const, + const struct unit_system *us, + struct swift_params *params, + const struct hydro_props *p, + const struct cosmology *cosmo) { + + /* Kernel properties */ + sp->eta_neighbours = parser_get_opt_param_float( + params, "Stars:resolution_eta", p->eta_neighbours); + + /* Tolerance for the smoothing length Newton-Raphson scheme */ + sp->h_tolerance = + parser_get_opt_param_float(params, "Stars:h_tolerance", p->h_tolerance); + + /* Get derived properties */ + sp->target_neighbours = pow_dimension(sp->eta_neighbours) * kernel_norm; + const float delta_eta = sp->eta_neighbours * (1.f + sp->h_tolerance); + sp->delta_neighbours = + (pow_dimension(delta_eta) - pow_dimension(sp->eta_neighbours)) * + kernel_norm; + + /* Number of iterations to converge h */ + sp->max_smoothing_iterations = parser_get_opt_param_int( + params, "Stars:max_ghost_iterations", p->max_smoothing_iterations); + + /* Time integration properties */ + const float max_volume_change = + parser_get_opt_param_float(params, "Stars:max_volume_change", -1); + if (max_volume_change == -1) + sp->log_max_h_change = p->log_max_h_change; + else + sp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv)); +} + +/** + * @brief Print the global properties of the stars scheme. + * + * @param sp The #stars_props. + */ +INLINE static void stars_props_print(const struct stars_props *sp) { + + /* Now stars */ + message("Stars kernel: %s with eta=%f (%.2f neighbours).", kernel_name, + sp->eta_neighbours, sp->target_neighbours); + + message("Stars relative tolerance in h: %.5f (+/- %.4f neighbours).", + sp->h_tolerance, sp->delta_neighbours); + + message( + "Stars integration: Max change of volume: %.2f " + "(max|dlog(h)/dt|=%f).", + pow_dimension(expf(sp->log_max_h_change)), sp->log_max_h_change); + + message("Maximal iterations in ghost task set to %d", + sp->max_smoothing_iterations); +} + +#if defined(HAVE_HDF5) +INLINE static void stars_props_print_snapshot(hid_t h_grpstars, + const struct stars_props *sp) { + + io_write_attribute_s(h_grpstars, "Kernel function", kernel_name); + io_write_attribute_f(h_grpstars, "Kernel target N_ngb", + sp->target_neighbours); + io_write_attribute_f(h_grpstars, "Kernel delta N_ngb", sp->delta_neighbours); + io_write_attribute_f(h_grpstars, "Kernel eta", sp->eta_neighbours); + io_write_attribute_f(h_grpstars, "Smoothing length tolerance", + sp->h_tolerance); + io_write_attribute_f(h_grpstars, "Volume log(max(delta h))", + sp->log_max_h_change); + io_write_attribute_f(h_grpstars, "Volume max change time-step", + pow_dimension(expf(sp->log_max_h_change))); + io_write_attribute_i(h_grpstars, "Max ghost iterations", + sp->max_smoothing_iterations); +} +#endif + +/** + * @brief Write a #stars_props struct to the given FILE as a stream of bytes. + * + * @param p the struct + * @param stream the file stream + */ +INLINE static void stars_props_struct_dump(const struct stars_props *p, + FILE *stream) { + restart_write_blocks((void *)p, sizeof(struct stars_props), 1, stream, + "starsprops", "stars props"); +} + +/** + * @brief Restore a stars_props struct from the given FILE as a stream of + * bytes. + * + * @param p the struct + * @param stream the file stream + */ +INLINE static void stars_props_struct_restore(const struct stars_props *p, + FILE *stream) { + restart_read_blocks((void *)p, sizeof(struct stars_props), 1, stream, NULL, + "stars props"); +} + +#endif /* SWIFT_GEAR_STAR_IO_H */ diff --git a/src/stars/GEAR/stars_part.h b/src/stars/GEAR/stars_part.h new file mode 100644 index 0000000000000000000000000000000000000000..bf68a580ef009f5a814fa17123301b5585d6084c --- /dev/null +++ b/src/stars/GEAR/stars_part.h @@ -0,0 +1,154 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Coypright (c) 2019 Loic Hausammann (loic.hausammann@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 + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_GEAR_STAR_PART_H +#define SWIFT_GEAR_STAR_PART_H + +/* Some standard headers. */ +#include <stdlib.h> + +/* Read additional subgrid models */ +#include "chemistry_struct.h" +#include "feedback_struct.h" +#include "tracers_struct.h" + +/** + * @brief Particle fields for the star particles. + * + * All quantities related to gravity are stored in the associate #gpart. + */ +struct spart { + + /*! Particle ID. */ + long long id; + + /*! Pointer to corresponding gravity part. */ + struct gpart* gpart; + + /*! Particle position. */ + double x[3]; + + /* Offset between current position and position at last tree rebuild. */ + float x_diff[3]; + + /* Offset between current position and position at last tree rebuild. */ + float x_diff_sort[3]; + + /*! Particle velocity. */ + float v[3]; + + /*! Star mass */ + float mass; + + /* Particle cutoff radius. */ + float h; + + /*! Union for the birth time and birth scale factor */ + union { + + /*! Birth time */ + float birth_time; + + /*! Birth scale factor */ + float birth_scale_factor; + }; + + /*! Particle time bin */ + timebin_t time_bin; + + struct { + + /* Number of neighbours. */ + float wcount; + + /* Number of neighbours spatial derivative. */ + float wcount_dh; + + } density; + + struct { + /*! birth density*/ + float density; + + /*! birth temperature*/ + float temperature; + + /*! birth mass */ + float mass; + } birth; + + /*! Feedback structure */ + struct feedback_spart_data feedback_data; + + /*! Tracer structure */ + struct tracers_xpart_data tracers_data; + + /*! Chemistry structure */ + struct chemistry_part_data chemistry_data; + +#ifdef SWIFT_DEBUG_CHECKS + + /* Time of the last drift */ + integertime_t ti_drift; + + /* Time of the last kick */ + integertime_t ti_kick; + +#endif + +#ifdef DEBUG_INTERACTIONS_STARS + /*! Number of interactions in the density SELF and PAIR */ + int num_ngb_density; + + /*! List of interacting particles in the density SELF and PAIR */ + long long ids_ngbs_density[MAX_NUM_OF_NEIGHBOURS_STARS]; + + /*! Number of interactions in the force SELF and PAIR */ + int num_ngb_force; + + /*! List of interacting particles in the force SELF and PAIR */ + long long ids_ngbs_force[MAX_NUM_OF_NEIGHBOURS_STARS]; +#endif + +} SWIFT_STRUCT_ALIGN; + +/** + * @brief Contains all the constants and parameters of the stars scheme + */ +struct stars_props { + + /*! Resolution parameter */ + float eta_neighbours; + + /*! Target weightd number of neighbours (for info only)*/ + float target_neighbours; + + /*! Smoothing length tolerance */ + float h_tolerance; + + /*! Tolerance on neighbour number (for info only)*/ + float delta_neighbours; + + /*! Maximal number of iterations to converge h */ + int max_smoothing_iterations; + + /*! Maximal change of h over one time-step */ + float log_max_h_change; +}; + +#endif /* SWIFT_GEAR_STAR_PART_H */