From e4e27f9b83b2bfa76755e2f41a065fd81f6ad3fa Mon Sep 17 00:00:00 2001
From: Loic Hausammann <loic_hausammann@hotmail.com>
Date: Fri, 13 Mar 2020 08:43:29 +0100
Subject: [PATCH] GEAR: fix the stellar feedback
---
src/feedback/GEAR/initial_mass_function.c | 2 +-
src/feedback/GEAR/interpolation.h | 10 ++++
src/feedback/GEAR/stellar_evolution.c | 72 +++++++++--------------
src/feedback/GEAR/stellar_evolution.h | 2 +-
src/feedback/GEAR/supernovae_ii.c | 3 +-
5 files changed, 41 insertions(+), 48 deletions(-)
diff --git a/src/feedback/GEAR/initial_mass_function.c b/src/feedback/GEAR/initial_mass_function.c
index 880bb73dbb..c0c5e2c81d 100644
--- a/src/feedback/GEAR/initial_mass_function.c
+++ b/src/feedback/GEAR/initial_mass_function.c
@@ -107,7 +107,7 @@ void initial_mass_function_integrate(const struct initial_mass_function *imf,
}
/* Integrate the data */
- while (m < imf->mass_limits[i + 1] && j < count) {
+ while ((m < imf->mass_limits[i + 1] || i == imf->n_parts - 1) && j < count) {
/* Compute the masses */
const float log_m1 = log_mass_min + (j - 1) * step_size;
diff --git a/src/feedback/GEAR/interpolation.h b/src/feedback/GEAR/interpolation.h
index 37706cda73..ed81e83898 100644
--- a/src/feedback/GEAR/interpolation.h
+++ b/src/feedback/GEAR/interpolation.h
@@ -28,6 +28,9 @@ enum interpolate_boundary_condition {
/* Zero (left boundary) and constant (right boundary) boundary conditions */
boundary_condition_zero_const,
+
+ /* constant boundary conditions */
+ boundary_condition_const,
};
struct interpolation_1d {
@@ -96,6 +99,9 @@ __attribute__((always_inline)) static INLINE void interpolate_1d_init(
case boundary_condition_zero_const:
interp->data[i] = 0;
break;
+ case boundary_condition_const:
+ interp->data[i] = data[0];
+ break;
default:
error("Interpolation type not implemented");
}
@@ -109,6 +115,7 @@ __attribute__((always_inline)) static INLINE void interpolate_1d_init(
interp->data[i] = 0;
break;
case boundary_condition_zero_const:
+ case boundary_condition_const:
interp->data[i] = interp->data[i - 1];
break;
default:
@@ -149,6 +156,8 @@ __attribute__((always_inline)) static INLINE float interpolate_1d(
case boundary_condition_zero:
case boundary_condition_zero_const:
return 0;
+ case boundary_condition_const:
+ return interp->data[0];
default:
error("Interpolation type not implemented");
}
@@ -160,6 +169,7 @@ __attribute__((always_inline)) static INLINE float interpolate_1d(
case boundary_condition_zero:
return 0;
case boundary_condition_zero_const:
+ case boundary_condition_const:
return interp->data[interp->N - 1];
default:
error("Interpolation type not implemented");
diff --git a/src/feedback/GEAR/stellar_evolution.c b/src/feedback/GEAR/stellar_evolution.c
index c4017f6ce8..c21d8f7e87 100644
--- a/src/feedback/GEAR/stellar_evolution.c
+++ b/src/feedback/GEAR/stellar_evolution.c
@@ -94,27 +94,28 @@ int stellar_evolution_compute_integer_number_supernovae(
* (solMass)
* @param m_end_step Mass of a star ending its life at the end of the step
* (solMass)
- * @param number_snia Number of SNIa produced by the stellar particle.
- * @param number_snii Number of SNII produced by the stellar particle.
+ * @param number_snia_f (Floating) Number of SNIa produced by the stellar particle.
+ * @param number_snii_f (Floating) Number of SNII produced by the stellar particle.
*
*/
void stellar_evolution_compute_continuous_feedback_properties(
struct spart* restrict sp, const struct stellar_model* sm,
const struct phys_const* phys_const, const float log_m_beg_step,
const float log_m_end_step, const float m_beg_step, const float m_end_step,
- const float m_init, const int number_snia, const int number_snii) {
+ const float m_init, const float number_snia_f, const float number_snii_f) {
/* Compute the mass ejected */
/* SNIa */
const float mass_snia =
- supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia;
+ supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia_f;
/* SNII */
const float mass_frac_snii = supernovae_ii_get_ejected_mass_fraction_processed(
&sm->snii, log_m_end_step, log_m_beg_step);
/* Sum the contributions from SNIa and SNII */
- sp->feedback_data.mass_ejected = mass_frac_snii * sp->birth.mass + mass_snia;
+ sp->feedback_data.mass_ejected = mass_frac_snii * sp->birth.mass
+ + mass_snia * phys_const->const_solar_mass;
if (sp->mass <= sp->feedback_data.mass_ejected) {
error("Stars cannot have negative mass. (%g <= %g). Initial mass = %g",
@@ -154,7 +155,7 @@ void stellar_evolution_compute_continuous_feedback_properties(
/* Add the Supernovae Ia */
sp->feedback_data.metal_mass_ejected[i] +=
- snia_yields[i] * number_snia * phys_const->const_solar_mass;
+ snia_yields[i] * number_snia_f * phys_const->const_solar_mass;
}
}
@@ -185,18 +186,15 @@ void stellar_evolution_compute_discrete_feedback_properties(
/* Get the normalization to the average */
const float normalization =
- number_snii == 0
- ? 0.
- : number_snii /
- (supernovae_ii_get_number_per_unit_mass(&sm->snii, m_end_step, m_beg_step) *
- m_init);
+ number_snii == 0 ? 0. :
+ number_snii / (supernovae_ii_get_number_per_unit_mass(&sm->snii, m_end_step, m_beg_step) *
+ m_init);
/* Compute the mass ejected */
/* SNIa */
const float mass_snia =
- (number_snia == 0)
- ? 0
- : (supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia);
+ (number_snia == 0) ? 0 :
+ (supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia);
/* SNII */
const float mass_snii =
@@ -305,50 +303,36 @@ void stellar_evolution_evolve_spart(
const float m_init = sp->sf_data.birth_mass / phys_const->const_solar_mass;
/* Compute number of SNIa */
- int number_snia = 0;
- if (can_produce_snia) {
- /* Compute rates */
- const float number_snia_f =
- supernovae_ia_get_number_per_unit_mass(&sm->snia, m_end_step, m_beg_step) * m_init;
-
- /* Get the integer number of supernovae */
- number_snia = stellar_evolution_compute_integer_number_supernovae(
- sp, number_snia_f, ti_begin, random_number_stellar_feedback_1);
-
- if (number_snia != 0) {
- message("Exploding %i SNIa", number_snia);
- }
- }
+ const float number_snia_f = can_produce_snia ?
+ supernovae_ia_get_number_per_unit_mass(&sm->snia, m_end_step, m_beg_step) * m_init : 0;
/* Compute number of SNII */
- int number_snii = 0;
- if (can_produce_snii) {
- /* Compute rates */
- const float number_snii_f =
- supernovae_ii_get_number_per_unit_mass(&sm->snii, m_end_step, m_beg_step) * m_init;
-
- /* Get the integer number of supernovae */
- number_snii = stellar_evolution_compute_integer_number_supernovae(
- sp, number_snii_f, ti_begin, random_number_stellar_feedback_2);
- if (number_snii != 0) {
- message("Exploding %i SNII", number_snii);
- }
- }
+ const float number_snii_f = can_produce_snii ?
+ supernovae_ii_get_number_per_unit_mass(&sm->snii, m_end_step, m_beg_step) * m_init : 0;
/* Does this star produce a supernovae? */
- if (number_snia == 0 && number_snii == 0) return;
+ if (number_snia_f == 0 && number_snii_f == 0) return;
- sp->feedback_data.number_sn = number_snia + number_snii;
+ sp->feedback_data.number_sn = number_snia_f + number_snii_f;
/* Compute the properties of the feedback (e.g. yields) */
if (sm->discrete_yields) {
+ /* Get the integer number of supernovae */
+ const int number_snia = stellar_evolution_compute_integer_number_supernovae(
+ sp, number_snia_f, ti_begin, random_number_stellar_feedback_1);
+
+ /* Get the integer number of supernovae */
+ const int number_snii = stellar_evolution_compute_integer_number_supernovae(
+ sp, number_snii_f, ti_begin, random_number_stellar_feedback_2);
+
stellar_evolution_compute_discrete_feedback_properties(
sp, sm, phys_const, log_m_beg_step, log_m_end_step, m_beg_step,
m_end_step, m_init, number_snia, number_snii);
+
} else {
stellar_evolution_compute_continuous_feedback_properties(
sp, sm, phys_const, log_m_beg_step, log_m_end_step, m_beg_step,
- m_end_step, m_init, number_snia, number_snii);
+ m_end_step, m_init, number_snia_f, number_snii_f);
}
}
diff --git a/src/feedback/GEAR/stellar_evolution.h b/src/feedback/GEAR/stellar_evolution.h
index 9ffda5ff67..82bbea7dde 100644
--- a/src/feedback/GEAR/stellar_evolution.h
+++ b/src/feedback/GEAR/stellar_evolution.h
@@ -39,7 +39,7 @@ void stellar_evolution_compute_continuous_feedback_properties(
struct spart* restrict sp, const struct stellar_model* sm,
const struct phys_const* phys_const, const float log_m_beg_step,
const float log_m_end_step, const float m_beg_step, const float m_end_step,
- const float m_init, const int number_snia, const int number_snii);
+ const float m_init, const float number_snia_f, const float number_snii_f);
void stellar_evolution_compute_discrete_feedback_properties(
struct spart* restrict sp, const struct stellar_model* sm,
const struct phys_const* phys_const, const float log_m_beg_step,
diff --git a/src/feedback/GEAR/supernovae_ii.c b/src/feedback/GEAR/supernovae_ii.c
index dc2a784c60..6220e97553 100644
--- a/src/feedback/GEAR/supernovae_ii.c
+++ b/src/feedback/GEAR/supernovae_ii.c
@@ -191,14 +191,13 @@ void supernovae_ii_read_yields_array(
/* Read the dataset */
io_read_array_dataset(group_id, hdf5_dataset_name, FLOAT, data, count);
- /* Integrate the yields */
initial_mass_function_integrate(&sm->imf, data, count, log_mass_min, step_size);
// TODO: decrease count in order to keep the same distance between points
/* Initialize the interpolation */
interpolate_1d_init(interp, log10(snii->mass_min), log10(snii->mass_max),
interpolation_size, log_mass_min, step_size, count, data,
- boundary_condition_zero_const);
+ boundary_condition_const);
/* Cleanup the memory */
--
GitLab