-
Matthieu Schaller authoredMatthieu Schaller authored
feedback_iact.h 11.89 KiB
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (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 <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_EAGLE_FEEDBACK_IACT_H
#define SWIFT_EAGLE_FEEDBACK_IACT_H
/* Local includes */
#include "random.h"
#include "timestep_sync_part.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 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;
/* 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;
}
/**
* @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 cosmo The cosmological model.
* @param ti_current Current integer time used value for seeding random number
* generator
*/
__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 struct cosmology *cosmo,
const integertime_t ti_current) {
#ifdef SWIFT_DEBUG_CHECKS
if (si->count_since_last_enrichment != 0)
error("Computing feedback from a star that should not");
#endif
/* 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);
/* 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
if (Omega_frac < 0. || Omega_frac > 1.01)
error(
"Invalid fraction of material to distribute for star ID=%lld "
"Omega_frac=%e count since last enrich=%d",
si->id, Omega_frac, si->count_since_last_enrichment);
#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;
const 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;
/* 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 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;
/* Update mass from SNIa */
const double delta_mass_from_SNIa =
si->feedback_data.to_distribute.mass_from_SNIa * Omega_frac;
pj->chemistry_data.mass_from_SNIa += delta_mass_from_SNIa;
/* 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 delta_mass_from_SNII =
si->feedback_data.to_distribute.mass_from_SNII * Omega_frac;
pj->chemistry_data.mass_from_SNII += delta_mass_from_SNII;
/* 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 delta_mass_from_AGB =
si->feedback_data.to_distribute.mass_from_AGB * Omega_frac;
pj->chemistry_data.mass_from_AGB += delta_mass_from_AGB;
/* 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;
/* 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
* (thermal SNIa + kinetic energy of ejecta + kinetic energy of star) */
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. */
const double new_thermal_energy = current_kinetic_energy_gas +
current_thermal_energy + injected_energy -
new_kinetic_energy_gas;
/* 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, u_new_enrich);
/* Finally, SNII stochastic feedback */
/* Get the SNII feedback properties */
const float prob = si->feedback_data.to_distribute.SNII_heating_probability;
/* Are we doing some SNII (big boys) feedback? */
if (prob > 0.f) {
/* Draw a random number (Note mixing both IDs) */
const float rand = random_unit_interval(si->id + pj->id, ti_current,
random_number_stellar_feedback);
/* Are we lucky? */
if (rand < prob) {
/* Compute new energy of this particle */
const double u_init = hydro_get_physical_internal_energy(pj, xpj, cosmo);
const float delta_u = si->feedback_data.to_distribute.SNII_delta_u;
const double u_new = u_init + delta_u;
/* Inject energy into the particle */
hydro_set_physical_internal_energy(pj, xpj, cosmo, u_new);
hydro_set_drifted_physical_internal_energy(pj, cosmo, u_new);
/* Impose maximal viscosity */
hydro_diffusive_feedback_reset(pj);
/* 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);
}
}
}
#endif /* SWIFT_EAGLE_FEEDBACK_IACT_H */