/******************************************************************************* * This file is part of SWIFT. * Copyright (c) 2018 Matthieu Schaller (schaller@strw.leidenuniv.nl) * * 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 . * ******************************************************************************/ #ifndef SWIFT_COLIBRE_FEEDBACK_IACT_H #define SWIFT_COLIBRE_FEEDBACK_IACT_H /* Local includes */ #include "event_logger.h" #include "feedback.h" #include "random.h" #include "rays.h" #include "timestep_sync_part.h" #include "tracers.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 xpj Extra particle data (not updated). * @param cosmo The cosmological model. * @param ti_current Current integer time value */ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_feedback_density(const float r2, const float *dx, const float hi, const float hj, struct spart *si, const struct part *pj, const struct xpart *xpj, const struct cosmology *cosmo, const struct feedback_props *fb_props, const integertime_t ti_current) { /* Get the gas mass. */ const float mj = hydro_get_mass(pj); /* Get r. */ const float r = sqrtf(r2); /* Compute the kernel function */ const float hi_inv = 1.0f / hi; const float ui = r * hi_inv; float wi; kernel_eval(ui, &wi); /* Add mass of pj to neighbour mass of si */ si->feedback_data.to_collect.ngb_mass += mj; si->feedback_data.to_collect.ngb_N++; /* Contribution to the star's surrounding gas density */ si->feedback_data.to_collect.ngb_rho += mj * wi; /* Neighbour internal energy */ const float uj = hydro_get_drifted_comoving_internal_energy(pj); /* Contribution to the star's surrounding internal energy */ si->feedback_data.to_collect.ngb_u += mj * uj * wi; const float Zj = chemistry_get_total_metal_mass_fraction_for_cooling(pj); /* Contribution to the star's surrounding metallicity (metal mass fraction */ si->feedback_data.to_collect.ngb_Z += mj * Zj * wi; /* Add contribution of pj to normalisation of density weighted fraction * which determines how much mass to distribute to neighbouring * gas particles */ const float rho = hydro_get_comoving_density(pj); if (rho != 0.f) si->feedback_data.to_collect.enrichment_weight_inv += wi / rho; /* Compute arc lengths in stellar isotropic feedback and collect * relevant data for later use in the feedback_apply loop */ /* Loop over rays */ for (int i = 0; i < colibre_SNII_feedback_num_of_rays; i++) { /* We generate two random numbers that we use * to randomly select the direction of the ith ray * and the associated mirror ray in SNII feedback */ /* Two random numbers in [0, 1[ */ const double rand_theta_SNII = random_unit_interval_part_ID_and_index( si->id, i, ti_current, random_number_isotropic_SNII_feedback_ray_theta); const double rand_phi_SNII = random_unit_interval_part_ID_and_index( si->id, i, ti_current, random_number_isotropic_SNII_feedback_ray_phi); /* Compute arclength for the true particle (SNII feedback) */ ray_minimise_arclength(dx, r, si->feedback_data.SNII_rays_true + i, ray_feedback_kinetic_true, pj->id, rand_theta_SNII, rand_phi_SNII, mj, si->feedback_data.SNII_rays_ext_true + i, pj->v); /* Compute arclength for the mirror particle (SNII feedback) */ ray_minimise_arclength(dx, r, si->feedback_data.SNII_rays_mirr + i, ray_feedback_kinetic_mirr, pj->id, rand_theta_SNII, rand_phi_SNII, mj, si->feedback_data.SNII_rays_ext_mirr + i, pj->v); } for (int i = 0; i < colibre_SNIa_feedback_num_of_rays; i++) { /* We generate two random numbers that we use * to randomly select the direction of the ith ray * in SNIa feedback */ /* Two random numbers in [0, 1[ */ const double rand_theta_SNIa = random_unit_interval_part_ID_and_index( si->id, i, ti_current, random_number_isotropic_SNIa_feedback_ray_theta); const double rand_phi_SNIa = random_unit_interval_part_ID_and_index( si->id, i, ti_current, random_number_isotropic_SNIa_feedback_ray_phi); /* Compute arclength (SNIa feedback, only thermal feedback) */ ray_minimise_arclength(dx, r, si->feedback_data.SNIa_rays + i, /*ray_type=*/ray_feedback_thermal, pj->id, rand_theta_SNIa, rand_phi_SNIa, mj, /*ray_ext=*/NULL, /*v=*/NULL); } } /** * @brief Prep1 interaction between two particles (non-symmetric). * Used for updating properties of gas particles neighbouring a star particle. * * @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 xpj Extra particle data (not updated). * @param cosmo The cosmological model. * @param ti_current Current integer time value */ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_feedback_prep1(const float r2, const float *dx, const float hi, const float hj, const struct spart *si, struct part *pj, const struct xpart *xpj, const struct cosmology *cosmo, const integertime_t ti_current) { /* Get the the number of SNII kinetic energy injections per stellar * particle at this time-step */ const int N_of_SNII_kinetic_events = si->feedback_data.to_distribute.SNII_num_of_kinetic_energy_inj; /* Should we even bother? */ if (N_of_SNII_kinetic_events > 0) { /* Energy per pair of two gas particles */ const double energy_per_pair = si->feedback_data.to_distribute.SNII_E_kinetic / N_of_SNII_kinetic_events; /* Loop over the SNII kick events. In each event, two gas * particles are kicked in exactly the opposite directions. */ for (int i = 0; i < N_of_SNII_kinetic_events; i++) { /* Find the particle that is closest to the ith ray OR the ith mirror ray */ if (pj->id == si->feedback_data.SNII_rays_true[i].id_min_length || pj->id == si->feedback_data.SNII_rays_mirr[i].id_min_length) { /* If this spart has the largest energy in its SNII kinetic energy * reservoir among all sparts that want to kick gas particle pj in this * time-step, then the gas particle will collect the id of this spart. */ if (pj->feedback_data.SNII_max_E_kin_pair_avail < energy_per_pair) { /* Update the stellar id and maximal available kinetic energy */ pj->feedback_data.SNII_stellar_id_kick_allowed = si->id; pj->feedback_data.SNII_max_E_kin_pair_avail = energy_per_pair; } /* If two stars have the save energies for the kicks, the one with the * larger id has the priority (very rare case) */ else if (pj->feedback_data.SNII_max_E_kin_pair_avail == energy_per_pair && pj->feedback_data.SNII_stellar_id_kick_allowed < si->id) { pj->feedback_data.SNII_stellar_id_kick_allowed = si->id; } } } } } /** * @brief Prep2 interaction between two particles (non-symmetric). * Used for updating properties of stellar particles. * * @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 xpj Extra particle data (not updated). * @param cosmo The cosmological model. * @param ti_current Current integer time value */ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_feedback_prep2(const float r2, const float *dx, const float hi, const float hj, struct spart *si, const struct part *pj, const struct xpart *xpj, const struct cosmology *cosmo, const integertime_t ti_current) { /* Get the the number of SNII kinetic energy injections per stellar * particle at this time-step */ const int N_of_SNII_kinetic_events = si->feedback_data.to_distribute.SNII_num_of_kinetic_energy_inj; /* Should we even bother? */ if (N_of_SNII_kinetic_events > 0) { /* Energy per pair of two gas particles */ const double energy_per_pair = si->feedback_data.to_distribute.SNII_E_kinetic / N_of_SNII_kinetic_events; /* Loop over the SNII kick events. In each event, two gas * particles are kicked in exactly the opposite directions. */ for (int i = 0; i < N_of_SNII_kinetic_events; i++) { /* Find the particle that is closest to the ith ray */ if (pj->id == si->feedback_data.SNII_rays_true[i].id_min_length) { /* Does this gas particle want to be kicked by this stellar particle * via ray i? If so, store this information in the ith ray extra struct */ if (pj->feedback_data.SNII_stellar_id_kick_allowed == si->id) { si->feedback_data.SNII_rays_ext_true[i].status = ray_feedback_kick_allowed; /* Take out energy from the reservoir if the other particle in the * pair is kicked too. Note that we can't do this step in the fourth * loop because that one only updates properties of gas particles. */ if (si->feedback_data.SNII_rays_ext_mirr[i].status == ray_feedback_kick_allowed) { si->feedback_data.SNII_E_kinetic_residual -= energy_per_pair; si->feedback_data.SNII_num_of_kinetic_energy_inj_residual--; } /* If we are using maximum_number_of_rays > 1, then for a given spart, * as soon as we have found the first ray that points at this gas * part, we stop. Otherwise, the same spart might kick the same gas * part twice in the same time-step (or even more times). */ break; } } /* Same as above but for the mirror ith ray */ else if (pj->id == si->feedback_data.SNII_rays_mirr[i].id_min_length) { if (pj->feedback_data.SNII_stellar_id_kick_allowed == si->id) { si->feedback_data.SNII_rays_ext_mirr[i].status = ray_feedback_kick_allowed; /* Take out energy from the reservoir if the other particle in the * pair is kicked too */ if (si->feedback_data.SNII_rays_ext_true[i].status == ray_feedback_kick_allowed) { si->feedback_data.SNII_E_kinetic_residual -= energy_per_pair; si->feedback_data.SNII_num_of_kinetic_energy_inj_residual--; } break; } } } } } /** * @brief Feedback interaction between two particles (non-symmetric). * Used for updating properties of gas particles neighbouring a star particle * * @param r2 Comoving square distance between the two particles. * @param dx Comoving vector separating both particles (si - pj). * @param hi Comoving smoothing-length of particle i. * @param hj Comoving smoothing-length of particle j. * @param si First (star) particle (not updated). * @param pj Second (gas) particle. * @param xpj Extra particle data. * @param with_cosmology Are we doing a cosmological run? * @param cosmo The cosmological model. * @param ti_current Current integer time used value for seeding random number * generator * @param time current physical time in the simulation * @param step current step counter */ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_feedback_apply(const float r2, const float *dx, const float hi, const float hj, const struct spart *si, struct part *pj, struct xpart *xpj, const int with_cosmology, const struct cosmology *cosmo, const struct hydro_props *hydro_props, const struct feedback_props *fb_props, const integertime_t ti_current, const double time, const int step) { /* Get r and 1/r. */ const float r = sqrtf(r2); const float r_inv = 1.f / r; /* Compute the kernel function */ const float hi_inv = 1.0f / hi; const float ui = r * hi_inv; float wi; kernel_eval(ui, &wi); /* Gas particle density */ const float rho_j = hydro_get_comoving_density(pj); /* Compute weighting for distributing feedback quantities */ float Omega_frac; if (rho_j != 0.f) { Omega_frac = si->feedback_data.to_distribute.enrichment_weight * wi / rho_j; } else { Omega_frac = 0.f; } // #ifdef SWIFT_DEBUG_CHECKS // MATTHIEU: Temporary check to track down possbile NaNs in dust. if (Omega_frac < 0. || Omega_frac > 1.01) warning( "Enrichment warning: star pid %lld -> gas pid %lld : Omega_frac=%e " "rho_j=%e rho_born=%e a=%e a_born=%e \n", si->id, pj->id, Omega_frac, rho_j, si->birth_density, cosmo->a, si->birth_time); // #endif /* Update particle mass */ const double current_mass = hydro_get_mass(pj); const double delta_mass = si->feedback_data.to_distribute.mass * Omega_frac; const double new_mass = current_mass + delta_mass; hydro_set_mass(pj, new_mass); /* Inverse of the new mass */ const double new_mass_inv = 1. / new_mass; /* Update total metallicity */ const double current_metal_mass_total = pj->chemistry_data.metal_mass_fraction_total * current_mass; const double delta_metal_mass_total = si->feedback_data.to_distribute.total_metal_mass * Omega_frac; double new_metal_mass_total = current_metal_mass_total + delta_metal_mass_total; pj->chemistry_data.metal_mass_fraction_total = new_metal_mass_total * new_mass_inv; /* Take ztime, it is either redshift or time, depending on with_cosmology */ double ztime; if (with_cosmology) { ztime = cosmo->z; } else { ztime = time; } /* Calculate mean metal weighted redshift */ if (pj->chemistry_data.metal_weighted_redshift < 0.f) pj->chemistry_data.metal_weighted_redshift = 0.f; if (delta_metal_mass_total > 0.f) { pj->chemistry_data.metal_weighted_redshift = ztime * delta_metal_mass_total + pj->chemistry_data.metal_weighted_redshift * pj->chemistry_data.track_of_metal_mass_total + pj->chemistry_data.metal_diffused_redshift; pj->chemistry_data.track_of_metal_mass_total += delta_metal_mass_total; pj->chemistry_data.metal_weighted_redshift /= pj->chemistry_data.track_of_metal_mass_total; /* Reset values */ pj->chemistry_data.metal_diffused_redshift = 0.f; } /* Update mass fraction of each tracked element */ for (int elem = 0; elem < chemistry_element_count; elem++) { const double current_metal_mass = pj->chemistry_data.metal_mass_fraction[elem] * current_mass; const double delta_metal_mass = si->feedback_data.to_distribute.metal_mass[elem] * Omega_frac; const double new_metal_mass = current_metal_mass + delta_metal_mass; pj->chemistry_data.metal_mass_fraction[elem] = new_metal_mass * new_mass_inv; } /* Update europium masses from neutron star mergers (NSM), common-envelop jets * SN (CEJSN) and collapsars */ const double current_mass_from_NSM = pj->chemistry_data.mass_fraction_from_NSM * current_mass; const double delta_mass_from_NSM = si->feedback_data.to_distribute.mass_from_NSM * Omega_frac; const double new_mass_from_NSM = current_mass_from_NSM + delta_mass_from_NSM; pj->chemistry_data.mass_fraction_from_NSM = new_mass_from_NSM * new_mass_inv; const double current_mass_from_CEJSN = pj->chemistry_data.mass_fraction_from_CEJSN * current_mass; const double delta_mass_from_CEJSN = si->feedback_data.to_distribute.mass_from_CEJSN * Omega_frac; const double new_mass_from_CEJSN = current_mass_from_CEJSN + delta_mass_from_CEJSN; pj->chemistry_data.mass_fraction_from_CEJSN = new_mass_from_CEJSN * new_mass_inv; const double current_mass_from_collapsar = pj->chemistry_data.mass_fraction_from_collapsar * current_mass; const double delta_mass_from_collapsar = si->feedback_data.to_distribute.mass_from_collapsar * Omega_frac; const double new_mass_from_collapsar = current_mass_from_collapsar + delta_mass_from_collapsar; pj->chemistry_data.mass_fraction_from_collapsar = new_mass_from_collapsar * new_mass_inv; #if !defined(DUST_NONE) /* Update mass fraction of each tracked dust grain species */ for (int grain = 0; grain < dust_grain_species_count; grain++) { const double current_grain_mass = pj->dust_data.grain_mass_fraction[grain] * current_mass; const double delta_grain_mass = si->feedback_data.to_distribute.dust_mass[grain] * Omega_frac; const double new_grain_mass = current_grain_mass + delta_grain_mass; pj->dust_data.grain_mass_fraction[grain] = new_grain_mass * new_mass_inv; } /* ! DO THIS BETTER -update mass fraction of depleted Fe from SNIa */ pj->dust_data.grain_iron_mass_fraction_from_SNIa = pj->dust_data.grain_iron_mass_fraction_from_SNIa * current_mass * new_mass_inv; #endif /* Update iron mass fraction from SNIa */ const double current_iron_from_SNIa_mass = pj->chemistry_data.iron_mass_fraction_from_SNIa * current_mass; const double delta_iron_from_SNIa_mass = si->feedback_data.to_distribute.Fe_mass_from_SNIa * Omega_frac; const double new_iron_from_SNIa_mass = current_iron_from_SNIa_mass + delta_iron_from_SNIa_mass; pj->chemistry_data.iron_mass_fraction_from_SNIa = new_iron_from_SNIa_mass * new_mass_inv; /* Calculate mean iron weighted redshift */ if (pj->chemistry_data.iron_weighted_redshift < 0.f) pj->chemistry_data.iron_weighted_redshift = 0.f; const double delta_iron_mass = si->feedback_data.to_distribute.metal_mass[chemistry_element_Fe] * Omega_frac; if (delta_iron_mass > 0.f) { pj->chemistry_data.iron_weighted_redshift = ztime * delta_iron_mass + pj->chemistry_data.iron_weighted_redshift * pj->chemistry_data.track_of_iron_mass + pj->chemistry_data.iron_diffused_redshift; pj->chemistry_data.track_of_iron_mass += delta_iron_mass; pj->chemistry_data.iron_weighted_redshift /= pj->chemistry_data.track_of_iron_mass; pj->chemistry_data.iron_diffused_redshift = 0.f; } /* Update mass from SNIa */ const double current_mass_from_SNIa = pj->chemistry_data.mass_fraction_from_SNIa * current_mass; const double delta_mass_from_SNIa = si->feedback_data.to_distribute.mass_from_SNIa * Omega_frac; const double new_mass_from_SNIa = current_mass_from_SNIa + delta_mass_from_SNIa; pj->chemistry_data.mass_fraction_from_SNIa = new_mass_from_SNIa * new_mass_inv; /* Update metal mass fraction from SNIa */ const double current_metal_mass_from_SNIa = pj->chemistry_data.metal_mass_fraction_from_SNIa * current_mass; const double delta_metal_mass_from_SNIa = si->feedback_data.to_distribute.metal_mass_from_SNIa * Omega_frac; const double new_metal_mass_from_SNIa = current_metal_mass_from_SNIa + delta_metal_mass_from_SNIa; pj->chemistry_data.metal_mass_fraction_from_SNIa = new_metal_mass_from_SNIa * new_mass_inv; /* Update mass from SNII */ const double current_mass_from_SNII = pj->chemistry_data.mass_fraction_from_SNII * current_mass; const double delta_mass_from_SNII = si->feedback_data.to_distribute.mass_from_SNII * Omega_frac; const double new_mass_from_SNII = current_mass_from_SNII + delta_mass_from_SNII; pj->chemistry_data.mass_fraction_from_SNII = new_mass_from_SNII * new_mass_inv; /* Update metal mass fraction from SNII */ const double current_metal_mass_from_SNII = pj->chemistry_data.metal_mass_fraction_from_SNII * current_mass; const double delta_metal_mass_from_SNII = si->feedback_data.to_distribute.metal_mass_from_SNII * Omega_frac; const double new_metal_mass_from_SNII = current_metal_mass_from_SNII + delta_metal_mass_from_SNII; pj->chemistry_data.metal_mass_fraction_from_SNII = new_metal_mass_from_SNII * new_mass_inv; /* Update mass from AGB */ const double current_mass_from_AGB = pj->chemistry_data.mass_fraction_from_AGB * current_mass; const double delta_mass_from_AGB = si->feedback_data.to_distribute.mass_from_AGB * Omega_frac; const double new_mass_from_AGB = current_mass_from_AGB + delta_mass_from_AGB; pj->chemistry_data.mass_fraction_from_AGB = new_mass_from_AGB * new_mass_inv; /* Update metal mass fraction from AGB */ const double current_metal_mass_from_AGB = pj->chemistry_data.metal_mass_fraction_from_AGB * current_mass; const double delta_metal_mass_from_AGB = si->feedback_data.to_distribute.metal_mass_from_AGB * Omega_frac; const double new_metal_mass_from_AGB = current_metal_mass_from_AGB + delta_metal_mass_from_AGB; pj->chemistry_data.metal_mass_fraction_from_AGB = new_metal_mass_from_AGB * new_mass_inv; /* SNII stochastic kinetic feedback begins. * To conserve linear momentum, it is done before the particle velocity * is recomputed due to the change in particle mass */ /* Get the the number of SNII kinetic energy injections per stellar * particle at this time-step */ const int N_of_SNII_kinetic_events = si->feedback_data.to_distribute.SNII_num_of_kinetic_energy_inj; /* Are we doing some SNII (big boys) kinetic feedback? */ if (N_of_SNII_kinetic_events > 0) { /* Loop over the number of SNII kinetic events. In each event, two * particles are kicked in exactly the opposite directions. */ for (int i = 0; i < N_of_SNII_kinetic_events; i++) { /* Find the particle that is closest to the ith ray OR the ith mirror ray */ if (pj->id == si->feedback_data.SNII_rays_true[i].id_min_length || pj->id == si->feedback_data.SNII_rays_mirr[i].id_min_length) { /* Are we allowed to kick? */ if (si->feedback_data.SNII_rays_ext_true[i].status == ray_feedback_kick_allowed && si->feedback_data.SNII_rays_ext_mirr[i].status == ray_feedback_kick_allowed) { /* Which particles have we caught: the original or the mirror one? */ const ray_feedback_type ray_type = (pj->id == si->feedback_data.SNII_rays_mirr[i].id_min_length) ? ray_feedback_kinetic_mirr : ray_feedback_kinetic_true; /* Get the SNII feedback kinetic energy per pair. We thus divide the * total kinetic energy we have from the *si stellar particle by the * number of events (in each event, two particles are kicked) */ const double energy_per_pair = si->feedback_data.to_distribute.SNII_E_kinetic / N_of_SNII_kinetic_events; /* Two random numbers in [0, 1[ */ const double rand_theta = random_unit_interval_part_ID_and_index( si->id, i, ti_current, random_number_isotropic_SNII_feedback_ray_theta); const double rand_phi = random_unit_interval_part_ID_and_index( si->id, i, ti_current, random_number_isotropic_SNII_feedback_ray_phi); /* Initialise the kick velocity vector and its modulus */ float v_kick[3] = {0.f, 0.f, 0.f}; float v_kick_abs = 0.f; const double mass_true = si->feedback_data.SNII_rays_true[i].mass; const double mass_mirr = si->feedback_data.SNII_rays_mirr[i].mass; /* Mass = 0 means that the ray does not point to any gas particle. * We need to check it for both the original and mirror ray. * We also make sure the energy is positive to avoid the possibility * of division by zero in the function below */ if (mass_true > 0.0 && mass_mirr > 0.0 && energy_per_pair > 0.0) { /* Compute the physical kick velocity in internal units */ ray_kinetic_feedback_compute_kick_velocity( v_kick, &v_kick_abs, si->feedback_data.SNII_rays_ext_true + i, si->feedback_data.SNII_rays_ext_mirr + i, ray_type, energy_per_pair, cosmo, current_mass, si->v, rand_theta, rand_phi, mass_true, mass_mirr); } /* Do the kicks by updating the particle velocity. * Note that xpj->v_full = a^2 * dx/dt, with x the comoving * coordinate. Therefore, a physical kick, dv, gets translated into a * code velocity kick, a * dv */ xpj->v_full[0] += v_kick[0] * cosmo->a; xpj->v_full[1] += v_kick[1] * cosmo->a; xpj->v_full[2] += v_kick[2] * cosmo->a; /* Update the signal velocity of the particle based on the velocity * kick */ hydro_set_v_sig_based_on_velocity_kick(pj, cosmo, v_kick_abs); /* Save relevant info for diagnostics of SNII kinetic feedback */ tracers_after_SNII_kinetic_feedback(xpj, with_cosmology, cosmo->a, time, v_kick_abs); /* Synchronize the particle on the timeline */ timestep_sync_part(pj); } } } } /* SNII stochastic kinetic feedback ends */ /* Now account in the fully energy and momentum conserving way * for the change in gas particle mass, energy and momentum due to * ABG feedback energy and stellar ejecta (with the mass contributed at this * time-step by all available feedback channels) moving at the star's velocity */ /* Compute the current kinetic energy */ const double current_v2 = xpj->v_full[0] * xpj->v_full[0] + xpj->v_full[1] * xpj->v_full[1] + xpj->v_full[2] * xpj->v_full[2]; const double current_kinetic_energy_gas = 0.5 * cosmo->a2_inv * current_mass * current_v2; /* Compute the current thermal energy */ const double current_thermal_energy = current_mass * hydro_get_physical_internal_energy(pj, xpj, cosmo); /* Apply conservation of momentum */ /* Update velocity following change in gas mass */ xpj->v_full[0] *= current_mass * new_mass_inv; xpj->v_full[1] *= current_mass * new_mass_inv; xpj->v_full[2] *= current_mass * new_mass_inv; /* Update velocity following addition of mass with different momentum */ xpj->v_full[0] += delta_mass * new_mass_inv * si->v[0]; xpj->v_full[1] += delta_mass * new_mass_inv * si->v[1]; xpj->v_full[2] += delta_mass * new_mass_inv * si->v[2]; /* Compute the new kinetic energy */ const double new_v2 = xpj->v_full[0] * xpj->v_full[0] + xpj->v_full[1] * xpj->v_full[1] + xpj->v_full[2] * xpj->v_full[2]; const double new_kinetic_energy_gas = 0.5 * cosmo->a2_inv * new_mass * new_v2; /* Energy injected * (kinetic energy of ejecta + kinetic energy of star AGB) */ const double injected_energy = si->feedback_data.to_distribute.energy * Omega_frac; /* Apply energy conservation to recover the new thermal energy of the gas * Note: in some specific cases the new_thermal_energy could be lower * than the current_thermal_energy, this is mainly the case if the change * in mass is relatively small and the velocity vectors between both the * gas particle and the star particle have a small angle. */ double new_thermal_energy = current_kinetic_energy_gas + current_thermal_energy + injected_energy - new_kinetic_energy_gas; /* In rare configurations the new thermal energy could become negative. * We must prevent that even if that implies a slight violation of the * conservation of total energy. */ /* The minimum energy (in units of energy not energy per mass) is * the total particle mass (including the mass to distribute) at the * minimal internal energy per unit mass */ const double min_u = hydro_props->minimal_internal_energy * new_mass; new_thermal_energy = max(new_thermal_energy, min_u); /* Convert this to a specific thermal energy */ const double u_new_enrich = new_thermal_energy * new_mass_inv; /* Do the energy injection. */ hydro_set_physical_internal_energy(pj, xpj, cosmo, u_new_enrich); hydro_set_drifted_physical_internal_energy(pj, cosmo, /*pfloor=*/NULL, u_new_enrich); /* Do SNIa and SNII stochastic thermal feedback */ /* Number of energy injections received by the gas particle in SNII and SNIa stochastic isotropic thermal feedback models */ int N_of_SNII_energy_inj_received_by_gas = 0; int N_of_SNIa_energy_inj_received_by_gas = 0; /* SNII stochastic thermal feedback begins */ /* Get the total number of SNII thermal energy injections per stellar particle * at this time-step */ const int N_of_SNII_thermal_energy_inj = si->feedback_data.to_distribute.SNII_num_of_thermal_energy_inj; /* Are we doing some SNII (big boys) thermal feedback? */ if (N_of_SNII_thermal_energy_inj > 0) { /* Maximum number of rays used in stellar feedback */ const int num_rays = colibre_SNII_feedback_num_of_rays; /* Find out how many rays this gas particle has received. * Note that this loop goes in the opposite direction compared to that * in SNII kicks (i.e. i-- in place of i++). That's because if we have more * than one ray per stellar particle per time-step, we then want to avoid * the situation in which the gas particle that is kicked is also heated */ for (int i = num_rays; i > num_rays - N_of_SNII_thermal_energy_inj; i--) { if (pj->id == si->feedback_data.SNII_rays_true[i - 1].id_min_length) N_of_SNII_energy_inj_received_by_gas++; } /* If the number of SNII energy injections > 0, do SNII feedback */ if (N_of_SNII_energy_inj_received_by_gas > 0) { /* Compute new energy of this particle */ const double u_init_SNII = hydro_get_physical_internal_energy(pj, xpj, cosmo); /* The total SNII energy the particle receives is the energy of one SNII * energy injection times the number of SNII energy injections. * N_of_SNII_energy_inj_received_by_gas is equal to the number of rays to * which the particle was found to be closest. * Since the heating probability << 1 for \Delta T=10^{7.5} K, * in practice N_of_SNII_energy_inj_received_by_gas is either 0 or 1 */ const float delta_u_SNII = si->feedback_data.to_distribute.SNII_delta_u * (float)N_of_SNII_energy_inj_received_by_gas; const double u_new_SNII = u_init_SNII + delta_u_SNII; #ifdef SWIFT_DEBUG_CHECKS message("SNII event at star age [Myr] = %.4f", si->feedback_data.to_distribute.SNII_star_age_Myr); #endif /* Inject energy into the particle */ hydro_set_physical_internal_energy(pj, xpj, cosmo, u_new_SNII); hydro_set_drifted_physical_internal_energy(pj, cosmo, /*pfloor=*/NULL, u_new_SNII); /* Impose maximal viscosity */ hydro_diffusive_feedback_reset(pj); /* Update cooling properties. */ cooling_update_feedback_particle(xpj); /* Mark this particle has having been heated by supernova feedback */ tracers_after_SNII_thermal_feedback( pj, xpj, with_cosmology, cosmo->a, time, delta_u_SNII / fb_props->temp_to_u_factor); /* Write the event to the SNII log file */ event_logger_SNII_log_event(si, pj, xpj, cosmo, si->SNII_f_E); #ifdef SWIFT_DEBUG_CHECKS event_logger_SNII_log_event_debug(time, si, pj, xpj, cosmo, step); #endif /* message( */ /* "We did some heating! id %llu star id %llu probability %.5e " */ /* "random_num %.5e du %.5e du/ini %.5e", */ /* pj->id, si->id, prob, rand, delta_u, delta_u / u_init); */ /* Synchronize the particle on the timeline */ timestep_sync_part(pj); } } /* SNIa stochastic thermal feedback begins */ /* Get the total number of SNIa thermal energy injections per stellar * particle at this time-step */ const int N_of_SNIa_thermal_energy_inj = si->feedback_data.to_distribute.SNIa_num_of_thermal_energy_inj; /* Are we doing some SNIa feedback? */ if (N_of_SNIa_thermal_energy_inj > 0) { /* Find out how many rays this gas particle has received. */ for (int i = 0; i < N_of_SNIa_thermal_energy_inj; i++) { if (pj->id == si->feedback_data.SNIa_rays[i].id_min_length) N_of_SNIa_energy_inj_received_by_gas++; } /* If the number of SNIa energy injections > 0, do SNIa feedback */ if (N_of_SNIa_energy_inj_received_by_gas > 0) { /* Compute new energy of this particle */ const double u_init_SNIa = hydro_get_physical_internal_energy(pj, xpj, cosmo); const float delta_u_SNIa = si->feedback_data.to_distribute.SNIa_delta_u; const double u_new_SNIa = u_init_SNIa + delta_u_SNIa * (float)N_of_SNIa_energy_inj_received_by_gas; /* Inject energy into the particle */ hydro_set_physical_internal_energy(pj, xpj, cosmo, u_new_SNIa); hydro_set_drifted_physical_internal_energy(pj, cosmo, /*pfloor=*/NULL, u_new_SNIa); /* Impose maximal viscosity */ hydro_diffusive_feedback_reset(pj); /* Update cooling properties. */ cooling_update_feedback_particle(xpj); /* Mark this particle has having been heated by supernova feedback */ tracers_after_SNIa_thermal_feedback( pj, xpj, with_cosmology, cosmo->a, time, delta_u_SNIa / fb_props->temp_to_u_factor); /* Write the event to the SNIa log file */ event_logger_SNIa_log_event(si, pj, xpj, cosmo); #ifdef SWIFT_DEBUG_CHECKS event_logger_SNIa_log_event_debug(time, si, pj, xpj, cosmo, step); #endif /* Synchronize the particle on the timeline */ timestep_sync_part(pj); } } /* We finish with the early feedback in the form of momentum injection. * We kick gas particle away from the star using the momentum available in the * timestep. This is done stochastically. * However, if delta_v is small enough (or even negative), this translates * into kicking all neighboring particles away from the star. */ const float delta_v = si->feedback_data.to_distribute.momentum_delta_v; const float momentum_prob = si->feedback_data.to_distribute.momentum_probability; /* Draw a random number (Note mixing both IDs) */ const float momentum_rand = random_unit_interval_two_IDs( si->id, pj->id, ti_current, random_number_stellar_winds); /* if lucky, perform the actual kick */ if (momentum_rand < momentum_prob) { /* Note that xpj->v_full = a^2 * dx/dt, with x the comoving coordinate. * Therefore, a physical kick, dv, gets translated into a * code velocity kick, a * dv */ xpj->v_full[0] -= delta_v * dx[0] * r_inv * cosmo->a; xpj->v_full[1] -= delta_v * dx[1] * r_inv * cosmo->a; xpj->v_full[2] -= delta_v * dx[2] * r_inv * cosmo->a; /* Store how much physical momentum is received, so we don't care about * cosmology */ xpj->tracers_data.stellar_wind_momentum_received += delta_v * current_mass; /* Mark this particle has having been heated by supernova feedback */ tracers_after_stellar_wind_feedback(xpj, with_cosmology, cosmo->a, time); /* Update the signal velocity of the particle based on the velocity kick */ hydro_set_v_sig_based_on_velocity_kick(pj, cosmo, delta_v); /* Synchronize the particle on the timeline */ timestep_sync_part(pj); } /* Put particles into HII regions */ const float HIIregion_prob = si->feedback_data.to_distribute.HIIregion_probability; /* Draw a random number (Note mixing both IDs) */ const float HIIregion_rand = random_unit_interval_two_IDs( si->id, pj->id, ti_current, random_number_HII_regions); /* if lucky, particle is now flagged as HII region */ if (HIIregion_rand < HIIregion_prob) { /* gas particle gets flagged as HII region */ xpj->tracers_data.HIIregion_timer_gas = si->feedback_data.to_distribute.HIIregion_endtime; xpj->tracers_data.HIIregion_starid = si->feedback_data.to_distribute.HIIregion_starid; /* If the particle hasn't received feedback from other sources, heat it up to the HII region energy level */ if (!N_of_SNII_energy_inj_received_by_gas && !N_of_SNIa_energy_inj_received_by_gas) { /* Compute new energy of this particle */ const double u_init_HII = hydro_get_physical_internal_energy(pj, xpj, cosmo); const float u_HII = si->feedback_data.to_distribute.HII_u; /* Put the particle on the HII temperature floor or leave it if it was already above */ if (u_init_HII < u_HII) { /* Inject energy into the particle */ hydro_set_physical_internal_energy(pj, xpj, cosmo, u_HII); hydro_set_drifted_physical_internal_energy(pj, cosmo, /*pfloor=*/NULL, u_HII); /* Make sure the particle does not cool any more */ hydro_set_physical_internal_energy_dt(pj, cosmo, 0.f); } /* Impose maximal viscosity */ hydro_diffusive_feedback_reset(pj); /* Synchronize the particle on the timeline */ timestep_sync_part(pj); } } } #endif /* SWIFT_COLIBRE_FEEDBACK_IACT_H */