/******************************************************************************* * This file is part of SWIFT. * Copyright (c) 2018 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 . * ******************************************************************************/ /* Include header */ #include "feedback.h" /* Local includes */ #include "cooling.h" #include "cosmology.h" #include "engine.h" #include "error.h" #include "feedback_properties.h" #include "hydro_properties.h" #include "part.h" #include "stellar_evolution.h" #include "units.h" #include /** * @brief Update the properties of the particle due to a supernovae. * * @param p The #part to consider. * @param xp The #xpart to consider. * @param e The #engine. */ void feedback_update_part(struct part* p, struct xpart* xp, const struct engine* e) { /* Did the particle receive a supernovae */ if (xp->feedback_data.delta_mass == 0) return; const struct cosmology* cosmo = e->cosmology; const struct pressure_floor_props* pressure_floor = e->pressure_floor_props; /* Turn off the cooling */ cooling_set_part_time_cooling_off(p, xp, e->time); /* Update mass */ const float old_mass = hydro_get_mass(p); const float new_mass = old_mass + xp->feedback_data.delta_mass; if (xp->feedback_data.delta_mass < 0.) { error("Delta mass smaller than 0"); } hydro_set_mass(p, new_mass); xp->feedback_data.delta_mass = 0; /* Update the density */ p->rho *= new_mass / old_mass; /* Update internal energy */ const float u = hydro_get_physical_internal_energy(p, xp, cosmo) * old_mass / new_mass; const float u_new = u + xp->feedback_data.delta_u; hydro_set_physical_internal_energy(p, xp, cosmo, u_new); hydro_set_drifted_physical_internal_energy(p, cosmo, pressure_floor, u_new); xp->feedback_data.delta_u = 0.; /* Update the velocities */ for (int i = 0; i < 3; i++) { const float dv = xp->feedback_data.delta_p[i] / new_mass; xp->v_full[i] += dv; p->v[i] += dv; xp->feedback_data.delta_p[i] = 0; } } /** * @brief Reset the gas particle-carried fields related to feedback at the * start of a step. * * Nothing to do here in the GEAR model. * * @param p The particle. * @param xp The extended data of the particle. */ void feedback_reset_part(struct part* p, struct xpart* xp) {} /** * @brief Compute the times for the stellar model. * * This function assumed to be called in the time step task. * * @param sp The #spart to act upon * @param with_cosmology Are we running with the cosmological expansion? * @param cosmo The current cosmological model. * @param star_age_beg_of_step (output) Age of the star at the beginning of the * step. * @param dt_enrichment (output) Time step for the stellar evolution. * @param ti_begin_star (output) Integer time at the beginning of the time step. * @param ti_current The current time (in integer) * @param time_base The time base. * @param time The current time (in double) */ void compute_time(struct spart* sp, const int with_cosmology, const struct cosmology* cosmo, double* star_age_beg_of_step, double* dt_enrichment, integertime_t* ti_begin_star, const integertime_t ti_current, const double time_base, const double time) { const integertime_t ti_step = get_integer_timestep(sp->time_bin); *ti_begin_star = get_integer_time_begin(ti_current, sp->time_bin); /* Get particle time-step */ double dt_star; if (with_cosmology) { dt_star = cosmology_get_delta_time(cosmo, *ti_begin_star, *ti_begin_star + ti_step); } else { dt_star = get_timestep(sp->time_bin, time_base); } /* Calculate age of the star at current time */ double star_age_end_of_step; if (with_cosmology) { if (cosmo->a > (double)sp->birth_scale_factor) star_age_end_of_step = cosmology_get_delta_time_from_scale_factors( cosmo, (double)sp->birth_scale_factor, cosmo->a); else star_age_end_of_step = 0.; } else { star_age_end_of_step = max(time - (double)sp->birth_time, 0.); } /* Get the length of the enrichment time-step */ *dt_enrichment = feedback_get_enrichment_timestep(sp, with_cosmology, cosmo, time, dt_star); *star_age_beg_of_step = star_age_end_of_step - *dt_enrichment; } /** * @brief Will this star particle want to do feedback during the next time-step? * * This is called in the time step task. * * In GEAR, we compute the full stellar evolution here. * * @param sp The particle to act upon * @param feedback_props The #feedback_props structure. * @param cosmo The current cosmological model. * @param us The unit system. * @param phys_const The #phys_const. * @param ti_current The current time (in integer) * @param time_base The time base. * @param time The physical time in internal units. */ void feedback_will_do_feedback( struct spart* sp, const struct feedback_props* feedback_props, const int with_cosmology, const struct cosmology* cosmo, const double time, const struct unit_system* us, const struct phys_const* phys_const, const integertime_t ti_current, const double time_base) { /* Zero the energy of supernovae */ sp->feedback_data.energy_ejected = 0; sp->feedback_data.will_do_feedback = 0; /* quit if the birth_scale_factor or birth_time is negative */ if (sp->birth_scale_factor < 0.0 || sp->birth_time < 0.0) return; /* Pick the correct table. (if only one table, threshold is < 0) */ const float metal = chemistry_get_star_total_iron_mass_fraction_for_feedback(sp); const float threshold = feedback_props->metallicity_max_first_stars; /* If metal < threshold, then sp is a first star particle. */ const int is_first_star = metal < threshold; const struct stellar_model* model = is_first_star ? &feedback_props->stellar_model_first_stars : &feedback_props->stellar_model; /* Compute the times */ double star_age_beg_step = 0; double dt_enrichment = 0; integertime_t ti_begin = 0; compute_time(sp, with_cosmology, cosmo, &star_age_beg_step, &dt_enrichment, &ti_begin, ti_current, time_base, time); #ifdef SWIFT_DEBUG_CHECKS if (sp->birth_time == -1.) error("Evolving a star particle that should not!"); if (star_age_beg_step + dt_enrichment < 0) { error("Negative age for a star"); } #endif /* Ensure that the age is positive (rounding errors) */ const double star_age_beg_step_safe = star_age_beg_step < 0 ? 0 : star_age_beg_step; /* A single star */ if (sp->star_type == single_star) { /* If the star has completely exploded, do not continue. This will also avoid NaN values in the liftetime if the mass is set to 0. Correction (28.04.2024): A bug fix in the mass of the star (see stellar_evolution.c in stellar_evolution_compute_X_feedback_properties, X=discrete, continuous) has changed the mass of the star from 0 to discrete_star_minimal_gravity_mass. Hence the fix is propagated here. */ if (sp->mass <= model->discrete_star_minimal_gravity_mass) { return; } /* Now, compute the stellar evolution state for individual star particles. */ stellar_evolution_evolve_individual_star(sp, model, cosmo, us, phys_const, ti_begin, star_age_beg_step_safe, dt_enrichment); } else { /* Compute the stellar evolution including SNe energy. This function treats the case of particles representing the whole IMF (star_type = star_population) and the particles representing only the continuous part of the IMF (star_type = star_population_continuous_IMF) */ stellar_evolution_evolve_spart(sp, model, cosmo, us, phys_const, ti_begin, star_age_beg_step_safe, dt_enrichment); } /* Apply the energy efficiency factor */ sp->feedback_data.energy_ejected *= feedback_props->supernovae_efficiency; /* Set the particle as doing some feedback */ sp->feedback_data.will_do_feedback = sp->feedback_data.energy_ejected != 0.; } /** * @brief Should this particle be doing any feedback-related operation? * * @param sp The #spart. * @param e The #engine. */ int feedback_is_active(const struct spart* sp, const struct engine* e) { /* the particle is inactive if its birth_scale_factor or birth_time is * negative */ if (sp->birth_scale_factor < 0.0 || sp->birth_time < 0.0) return 0; return sp->feedback_data.will_do_feedback; } /** * @brief Returns the length of time since the particle last did * enrichment/feedback. * * @param sp The #spart. * @param with_cosmology Are we running with cosmological time integration on? * @param cosmo The cosmological model. * @param time The current time (since the Big Bang / start of the run) in * internal units. * @param dt_star the length of this particle's time-step in internal units. * @return The length of the enrichment step in internal units. */ double feedback_get_enrichment_timestep(const struct spart* sp, const int with_cosmology, const struct cosmology* cosmo, const double time, const double dt_star) { return dt_star; } /** * @brief Prepares a s-particle for its feedback interactions * * @param sp The particle to act upon */ void feedback_init_spart(struct spart* sp) { sp->feedback_data.enrichment_weight = 0.f; } /** * @brief Prepare the feedback fields after a star is born. * * This function is called in the functions sink_copy_properties_to_star() and * star_formation_copy_properties(). * * @param sp The #spart to act upon. * @param feedback_props The feedback perties to use. * @param star_type The stellar particle type. */ void feedback_init_after_star_formation( struct spart* sp, const struct feedback_props* feedback_props, const enum stellar_type star_type) { feedback_init_spart(sp); /* Zero the energy of supernovae */ sp->feedback_data.energy_ejected = 0; /* Activate the feedback loop for the first step */ sp->feedback_data.will_do_feedback = 1; /* Give to the star its appropriate type: single star, continuous IMF star or single population star */ sp->star_type = star_type; } /** * @brief Prepares a star's feedback field before computing what * needs to be distributed. * * This is called in the stars ghost. */ void feedback_reset_feedback(struct spart* sp, const struct feedback_props* feedback_props) {} /** * @brief Initialises the s-particles feedback props 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 feedback_props The properties of the feedback model. */ void feedback_first_init_spart(struct spart* sp, const struct feedback_props* feedback_props) { feedback_init_spart(sp); /* Zero the energy of supernovae */ sp->feedback_data.energy_ejected = 0; /* Activate the feedback loop for the first step */ sp->feedback_data.will_do_feedback = 1; } /** * @brief Initialises the s-particles feedback props 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 feedback_props The properties of the feedback model. */ void feedback_prepare_spart(struct spart* sp, const struct feedback_props* feedback_props) {} /** * @brief Prepare a #spart for the feedback task. * * This is called in the stars ghost task. * * In here, we only need to add the missing coefficients. * * @param sp The particle to act upon * @param feedback_props The #feedback_props structure. * @param cosmo The current cosmological model. * @param us The unit system. * @param phys_const The #phys_const. * @param star_age_beg_step The age of the star at the star of the time-step in * internal units. * @param dt The time-step size of this star in internal units. * @param time The physical time in internal units. * @param ti_begin The integer time at the beginning of the step. * @param with_cosmology Are we running with cosmology on? */ void feedback_prepare_feedback(struct spart* restrict sp, const struct feedback_props* feedback_props, const struct cosmology* cosmo, const struct unit_system* us, const struct phys_const* phys_const, const double star_age_beg_step, const double dt, const double time, const integertime_t ti_begin, const int with_cosmology) { /* Add missing h factor */ const float hi_inv = 1.f / sp->h; const float hi_inv_dim = pow_dimension(hi_inv); /* 1/h^d */ sp->feedback_data.enrichment_weight *= hi_inv_dim; } /** * @brief Write a feedback struct to the given FILE as a stream of bytes. * * @param feedback the struct * @param stream the file stream */ void feedback_struct_dump(const struct feedback_props* feedback, FILE* stream) { restart_write_blocks((void*)feedback, sizeof(struct feedback_props), 1, stream, "feedback", "feedback function"); stellar_evolution_dump(&feedback->stellar_model, stream); if (feedback->metallicity_max_first_stars != -1) { stellar_evolution_dump(&feedback->stellar_model_first_stars, stream); } } /** * @brief Restore a feedback struct from the given FILE as a stream of * bytes. * * @param feedback the struct * @param stream the file stream */ void feedback_struct_restore(struct feedback_props* feedback, FILE* stream) { restart_read_blocks((void*)feedback, sizeof(struct feedback_props), 1, stream, NULL, "feedback function"); stellar_evolution_restore(&feedback->stellar_model, stream); if (feedback->metallicity_max_first_stars != -1) { stellar_evolution_restore(&feedback->stellar_model_first_stars, stream); } } /** * @brief Clean the allocated memory. * * @param feedback the #feedback_props. */ void feedback_clean(struct feedback_props* feedback) { stellar_evolution_clean(&feedback->stellar_model); if (feedback->metallicity_max_first_stars != -1) { stellar_evolution_clean(&feedback->stellar_model_first_stars); } }