diff --git a/configure.ac b/configure.ac
index 2a15c2190d68dfda1d517a98d77057d0741a669f..6dc5ffd28ee1a17197e42dce7c1b8451256025b5 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1329,7 +1329,7 @@ case "$with_subgrid" in
with_subgrid_entropy_floor=none
with_subgrid_stars=GEAR
with_subgrid_star_formation=GEAR
- with_subgrid_feedback=thermal
+ with_subgrid_feedback=none
;;
EAGLE)
with_subgrid_cooling=EAGLE
@@ -1338,7 +1338,7 @@ case "$with_subgrid" in
with_subgrid_entropy_floor=EAGLE
with_subgrid_stars=EAGLE
with_subgrid_star_formation=EAGLE
- with_subgrid_feedback=none
+ with_subgrid_feedback=EAGLE
;;
*)
AC_MSG_ERROR([Unknown subgrid choice: $with_subgrid])
@@ -1727,11 +1727,8 @@ case "$with_stars" in
GEAR)
AC_DEFINE([STARS_GEAR], [1], [GEAR stellar model])
;;
- EAGLE)
- AC_DEFINE([STARS_EAGLE], [1], [EAGLE stellar model])
- ;;
none)
- AC_DEFINE([STARS_NONE], [1], [None stellar model])
+ AC_DEFINE([STARS_NONE], [1], [Basic stellar model])
;;
*)
@@ -1742,7 +1739,7 @@ esac
# Feedback model
AC_ARG_WITH([feedback],
[AS_HELP_STRING([--with-feedback=<model>],
- [Feedback model to use @<:@none, thermal, debug default: none@:>@]
+ [Feedback model to use @<:@none, EAGLE, debug default: none@:>@]
)],
[with_feedback="$withval"],
[with_feedback="none"]
@@ -1757,10 +1754,11 @@ if test "$with_subgrid" != "none"; then
fi
case "$with_feedback" in
- thermal)
- AC_DEFINE([FEEDBACK_THERMAL], [1], [Thermal Blastwave])
+ EAGLE)
+ AC_DEFINE([FEEDBACK_EAGLE], [1], [EAGLE stellar feedback and evolution model])
;;
none)
+ AC_DEFINE([FEEDBACK_NONE], [1], [No feedback])
;;
*)
@@ -1889,6 +1887,9 @@ AM_CONDITIONAL([HAVE_DOXYGEN], [test "$ac_cv_path_ac_pt_DX_DOXYGEN" != ""])
# Check if using EAGLE cooling
AM_CONDITIONAL([HAVEEAGLECOOLING], [test $with_cooling = "EAGLE"])
+# Check if using EAGLE feedback
+AM_CONDITIONAL([HAVEEAGLEFEEDBACK], [test $with_feedback = "EAGLE"])
+
# Handle .in files.
AC_CONFIG_FILES([Makefile src/Makefile examples/Makefile examples/Cooling/CoolingRates/Makefile doc/Makefile doc/Doxyfile tests/Makefile])
AC_CONFIG_FILES([argparse/Makefile tools/Makefile])
diff --git a/src/Makefile.am b/src/Makefile.am
index d7e4249a7ff67132505e3a7df8a134d4cd8b266c..8509880970f9ccbef67ef141d63785590d7f51c5 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -59,6 +59,12 @@ if HAVEEAGLECOOLING
EAGLE_COOLING_SOURCES += cooling/EAGLE/cooling.c cooling/EAGLE/cooling_tables.c
endif
+# source files for EAGLE feedback
+EAGLE_FEEDBACK_SOURCES =
+if HAVEEAGLEFEEDBACK
+EAGLE_FEEDBACK_SOURCES += feedback/EAGLE/feedback.c
+endif
+
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c engine_maketasks.c \
engine_marktasks.c engine_drift.c serial_io.c timers.c debug.c scheduler.c \
@@ -71,7 +77,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c engine_maketasks.c
collectgroup.c hydro_space.c equation_of_state.c \
chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c \
outputlist.c velociraptor_dummy.c logger_io.c memuse.c \
- $(EAGLE_COOLING_SOURCES)
+ $(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES)
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
diff --git a/src/feedback.h b/src/feedback.h
new file mode 100644
index 0000000000000000000000000000000000000000..e4987f916822fd90affd27372a51b6f11cf73287
--- /dev/null
+++ b/src/feedback.h
@@ -0,0 +1,36 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_FEEDBACK_H
+#define SWIFT_FEEDBACK_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Select the correct feedback model */
+#if defined(FEEDBACK_NONE)
+#include "./feedback/Default/feedback.h"
+#include "./deedback/Default/feedback_iact.h"
+#elif defined(FEEDBACK_EAGLE)
+#include "./feedback/EAGLE/feedback.h"
+#include "./feedback/EAGLE/feedback_iact.h"
+#else
+#error "Invalid choice of feedback model"
+#endif
+
+#endif
diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c
new file mode 100644
index 0000000000000000000000000000000000000000..0b4f99ed5b690cbda8fc0b84922716506b1fe6ba
--- /dev/null
+++ b/src/feedback/EAGLE/feedback.c
@@ -0,0 +1,742 @@
+/*******************************************************************************
+ * 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 <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+
+#include "feedback.h"
+
+#include "hydro_properties.h"
+#include "imf.h"
+#include "interpolate.h"
+#include "yield_tables.h"
+
+/**
+ * @brief determine which metallicity bin star belongs to for AGB, compute bin
+ * indices and offsets
+ *
+ * @param iz_low Pointer to index of metallicity bin to which the star belongs
+ * (to be calculated in this function)
+ * @param iz_high Pointer to index of metallicity bin to above the star's
+ * metallicity (to be calculated in this function)
+ * @param dz metallicity bin offset
+ * @param log10_metallicity log10 of star metallicity
+ * @param star_properties stars_props data structure
+ */
+inline static void determine_bin_yield_AGB(
+ int* iz_low, int* iz_high, float* dz, float log10_metallicity,
+ const struct feedback_props* feedback_props) {
+
+ // MATTHIEU
+
+ /* if (log10_metallicity > log10_min_metallicity) { */
+ /* /\* Find metallicity bin which contains the star's metallicity *\/ */
+ /* int j; */
+ /* for (j = 0; j < star_properties->feedback.AGB_n_z - 1 && */
+ /* log10_metallicity > */
+ /* star_properties->feedback.yield_AGB.metallicity[j + 1]; */
+ /* j++) */
+ /* ; */
+ /* *iz_low = j; */
+ /* *iz_high = *iz_low + 1; */
+
+ /* /\* Compute offset *\/ */
+ /* if (log10_metallicity >= */
+ /* star_properties->feedback.yield_AGB.metallicity[0] && */
+ /* log10_metallicity <= */
+ /* star_properties->feedback.yield_AGB */
+ /* .metallicity[star_properties->feedback.AGB_n_z - 1]) */
+ /* *dz = log10_metallicity - */
+ /* star_properties->feedback.yield_AGB.metallicity[*iz_low]; */
+ /* else */
+ /* *dz = 0; */
+
+ /* /\* Normalize offset *\/ */
+ /* float deltaz = star_properties->feedback.yield_AGB.metallicity[*iz_high] - */
+ /* star_properties->feedback.yield_AGB.metallicity[*iz_low]; */
+
+ /* if (deltaz > 0) */
+ /* *dz /= deltaz; */
+ /* else */
+ /* dz = 0; */
+ /* } else { */
+ /* *iz_low = 0; */
+ /* *iz_high = 0; */
+ /* *dz = 0; */
+ /* } */
+}
+
+/**
+ * @brief determine which metallicity bin star belongs to for SNII, compute bin
+ * indices and offsets
+ *
+ * @param iz_low Pointer to index of metallicity bin to which the star belongs
+ * (to be calculated in this function)
+ * @param iz_high Pointer to index of metallicity bin to above the star's
+ * metallicity (to be calculated in this function)
+ * @param dz metallicity bin offset
+ * @param log10_metallicity log10 of star metallicity
+ * @param star_properties stars_props data structure
+ */
+inline static void determine_bin_yield_SNII(
+ int* iz_low, int* iz_high, float* dz, float log10_metallicity,
+ const struct feedback_props* feedback_props) {
+
+ // MATTHIEU
+
+ /* if (log10_metallicity > log10_min_metallicity) { */
+ /* /\* Find metallicity bin which contains the star's metallicity *\/ */
+ /* int j; */
+ /* for (j = 0; j < star_properties->feedback.SNII_n_z - 1 && */
+ /* log10_metallicity > */
+ /* star_properties->feedback.yield_SNII.metallicity[j + 1]; */
+ /* j++) */
+ /* ; */
+ /* *iz_low = j; */
+ /* *iz_high = *iz_low + 1; */
+
+ /* /\* Compute offset *\/ */
+ /* if (log10_metallicity >= */
+ /* star_properties->feedback.yield_SNII.metallicity[0] && */
+ /* log10_metallicity <= */
+ /* star_properties->feedback.yield_SNII */
+ /* .metallicity[star_properties->feedback.SNII_n_z - 1]) */
+ /* *dz = log10_metallicity - */
+ /* star_properties->feedback.yield_SNII.metallicity[*iz_low]; */
+ /* else */
+ /* *dz = 0; */
+
+ /* /\* Normalize offset *\/ */
+ /* float deltaz = star_properties->feedback.yield_SNII.metallicity[*iz_high] - */
+ /* star_properties->feedback.yield_SNII.metallicity[*iz_low]; */
+
+ /* if (deltaz > 0) */
+ /* *dz = *dz / deltaz; */
+ /* else */
+ /* dz = 0; */
+ /* } else { */
+ /* *iz_low = 0; */
+ /* *iz_high = 0; */
+ /* *dz = 0; */
+ /* } */
+}
+
+
+/**
+ * @brief compute enrichment and feedback due to SNIa. To do this compute the
+ * number of SNIa that occur during the timestep, multiply by constants read
+ * from tables.
+ *
+ * @param log10_min_mass log10 mass at the end of step
+ * @param log10_max_mass log10 mass at the beginning of step
+ * @param stars star properties data structure
+ * @param sp spart we are computing feedback from
+ * @param star_age_Gyr age of star in Gyr
+ * @param dt_Gyr timestep dt in Gyr
+ */
+inline static void evolve_SNIa(float log10_min_mass, float log10_max_mass,
+ const struct feedback_props* feedback_props,
+ struct spart* restrict sp, float star_age_Gyr,
+ float dt_Gyr) {
+
+ /* Check if we're outside the mass range for SNIa */
+ if (log10_min_mass >= feedback_props->log10_SNIa_max_mass_msun) return;
+
+ /* If the max mass is outside the mass range update it to be the maximum and
+ * use updated values for the star's age and timestep in this function */
+ if (log10_max_mass > feedback_props->log10_SNIa_max_mass_msun) {
+ log10_max_mass = feedback_props->log10_SNIa_max_mass_msun;
+ float lifetime_Gyr =
+ lifetime_in_Gyr(exp(M_LN10 * feedback_props->log10_SNIa_max_mass_msun),
+ sp->chemistry_data.metal_mass_fraction_total, feedback_props);
+ dt_Gyr = star_age_Gyr + dt_Gyr - lifetime_Gyr;
+ star_age_Gyr = lifetime_Gyr;
+ }
+
+ /* compute the number of SNIa */
+ /* Efolding (Forster 2006) */
+ float num_SNIa_per_msun =
+ feedback_props->SNIa_efficiency *
+ (exp(-star_age_Gyr / feedback_props->SNIa_timescale_Gyr) -
+ exp(-(star_age_Gyr + dt_Gyr) / feedback_props->SNIa_timescale_Gyr)) *
+ sp->mass_init;
+
+ sp->feedback_data.to_distribute.num_SNIa =
+ num_SNIa_per_msun / feedback_props->const_solar_mass;
+
+ /* compute mass fractions of each metal */
+ for (int i = 0; i < chemistry_element_count; i++) {
+ sp->feedback_data.to_distribute.metal_mass[i] +=
+ num_SNIa_per_msun * feedback_props->yield_SNIa_IMF_resampled[i];
+ }
+
+ /* Update the metallicity of the material released */
+ sp->feedback_data.to_distribute.metal_mass_from_SNIa +=
+ num_SNIa_per_msun * feedback_props->yield_SNIa_total_metals_IMF_resampled;
+
+ /* Update the metal mass produced */
+ sp->feedback_data.to_distribute.total_metal_mass +=
+ num_SNIa_per_msun * feedback_props->yield_SNIa_total_metals_IMF_resampled;
+
+ /* Compute the mass produced by SNIa */
+ sp->feedback_data.to_distribute.mass_from_SNIa +=
+ num_SNIa_per_msun * feedback_props->yield_SNIa_total_metals_IMF_resampled;
+
+ /* Compute the iron mass produced */
+ sp->feedback_data.to_distribute.Fe_mass_from_SNIa +=
+ num_SNIa_per_msun *
+ feedback_props->yield_SNIa_IMF_resampled[chemistry_element_Fe];
+}
+
+/**
+ * @brief compute enrichment and feedback due to SNII. To do this, integrate the
+ * IMF weighted by the yields read from tables for each of the quantities of
+ * interest.
+ *
+ * @param log10_min_mass log10 mass at the end of step
+ * @param log10_max_mass log10 mass at the beginning of step
+ * @param stellar_yields array to store calculated yields for passing to
+ * integrate_imf
+ * @param stars star properties data structure
+ * @param sp spart we are computing feedback from
+ */
+inline static void evolve_SNII(float log10_min_mass, float log10_max_mass,
+ float* stellar_yields,
+ const struct feedback_props* restrict feedback_props,
+ struct spart* restrict sp) {
+ int low_imf_mass_bin_index, high_imf_mass_bin_index, mass_bin_index;
+
+ /* If mass at beginning of step is less than tabulated lower bound for IMF,
+ * limit it.*/
+ if (log10_min_mass < feedback_props->log10_SNII_min_mass_msun)
+ log10_min_mass = feedback_props->log10_SNII_min_mass_msun;
+
+ /* If mass at end of step is greater than tabulated upper bound for IMF, limit
+ * it.*/
+ if (log10_max_mass > feedback_props->log10_SNII_max_mass_msun)
+ log10_max_mass = feedback_props->log10_SNII_max_mass_msun;
+
+ /* Don't do anything if the stellar mass hasn't decreased by the end of the
+ * step */
+ if (log10_min_mass >= log10_max_mass) return;
+
+ /* determine which IMF mass bins contribute to the integral */
+ determine_imf_bins(log10_min_mass, log10_max_mass, &low_imf_mass_bin_index,
+ &high_imf_mass_bin_index, feedback_props);
+
+ /* Integrate IMF to determine number of SNII */
+ sp->feedback_data.to_distribute.num_SNII =
+ integrate_imf(log10_min_mass, log10_max_mass, 0, stellar_yields, feedback_props);
+
+ /* determine which metallicity bin and offset this star belongs to */
+ int iz_low = 0, iz_high = 0, low_index_3d, high_index_3d, low_index_2d, high_index_2d;
+ float dz = 0.;
+ determine_bin_yield_SNII(&iz_low, &iz_high, &dz,
+ log10(sp->chemistry_data.metal_mass_fraction_total),
+ feedback_props);
+
+ /* compute metals produced */
+ float metal_mass_released[chemistry_element_count], metal_mass_released_total;
+ for (int elem = 0; elem < chemistry_element_count; elem++) {
+ for (mass_bin_index = low_imf_mass_bin_index;
+ mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) {
+ low_index_3d = feedback_row_major_index_3d(
+ iz_low, elem, mass_bin_index, feedback_props->SNII_n_z,
+ chemistry_element_count, feedback_props->n_imf_mass_bins);
+ high_index_3d = feedback_row_major_index_3d(
+ iz_high, elem, mass_bin_index, feedback_props->SNII_n_z,
+ chemistry_element_count, feedback_props->n_imf_mass_bins);
+ low_index_2d = feedback_row_major_index_2d(
+ iz_low, mass_bin_index, feedback_props->SNII_n_z,
+ feedback_props->n_imf_mass_bins);
+ high_index_2d = feedback_row_major_index_2d(
+ iz_high, mass_bin_index, feedback_props->SNII_n_z,
+ feedback_props->n_imf_mass_bins);
+ stellar_yields[mass_bin_index] =
+ (1 - dz) *
+ (feedback_props->yield_SNII.yield_IMF_resampled[low_index_3d] +
+ sp->chemistry_data.metal_mass_fraction[elem] *
+ feedback_props->yield_SNII
+ .ejecta_IMF_resampled[low_index_2d]) +
+ dz * (feedback_props->yield_SNII.yield_IMF_resampled[high_index_3d] +
+ sp->chemistry_data.metal_mass_fraction[elem] *
+ feedback_props->yield_SNII
+ .ejecta_IMF_resampled[high_index_2d]);
+ }
+
+ metal_mass_released[elem] =
+ integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, feedback_props);
+ }
+
+ /* Compute mass produced */
+ for (mass_bin_index = low_imf_mass_bin_index;
+ mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) {
+ low_index_2d = feedback_row_major_index_2d(iz_low, mass_bin_index,
+ feedback_props->SNII_n_z,
+ feedback_props->n_imf_mass_bins);
+ high_index_2d = feedback_row_major_index_2d(
+ iz_high, mass_bin_index, feedback_props->SNII_n_z,
+ feedback_props->n_imf_mass_bins);
+ stellar_yields[mass_bin_index] =
+ (1 - dz) * (feedback_props->yield_SNII
+ .total_metals_IMF_resampled[low_index_2d] +
+ sp->chemistry_data.metal_mass_fraction_total *
+ feedback_props->yield_SNII
+ .ejecta_IMF_resampled[low_index_2d]) +
+ dz * (feedback_props->yield_SNII
+ .total_metals_IMF_resampled[high_index_2d] +
+ sp->chemistry_data.metal_mass_fraction_total *
+ feedback_props->yield_SNII
+ .ejecta_IMF_resampled[high_index_2d]);
+ }
+
+ metal_mass_released_total =
+ integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, feedback_props);
+
+ /* yield normalization */
+ float mass_ejected, mass_released;
+
+ /* zero all negative values */
+ for (int i = 0; i < chemistry_element_count; i++)
+ metal_mass_released[i] = max(metal_mass_released[i], 0.f);
+
+ metal_mass_released_total = max(metal_mass_released_total, 0.f);
+
+ /* compute the total metal mass ejected from the star*/
+ for (mass_bin_index = low_imf_mass_bin_index;
+ mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) {
+ low_index_2d = feedback_row_major_index_2d(iz_low, mass_bin_index,
+ feedback_props->SNII_n_z,
+ feedback_props->n_imf_mass_bins);
+ high_index_2d = feedback_row_major_index_2d(
+ iz_high, mass_bin_index, feedback_props->SNII_n_z,
+ feedback_props->n_imf_mass_bins);
+ stellar_yields[mass_bin_index] =
+ (1 - dz) *
+ feedback_props->yield_SNII.ejecta_IMF_resampled[low_index_2d] +
+ dz * feedback_props->yield_SNII.ejecta_IMF_resampled[high_index_2d];
+ }
+
+ mass_ejected =
+ integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, feedback_props);
+
+ /* compute the total mass released */
+ mass_released = metal_mass_released_total +
+ metal_mass_released[chemistry_element_H] +
+ metal_mass_released[chemistry_element_He];
+
+ /* normalize the yields */
+ if (mass_released > 0) {
+ /* Set normalisation factor. Note additional multiplication by the star
+ * initial mass as tables are per initial mass */
+ const float norm_factor = sp->mass_init * mass_ejected / mass_released;
+
+ for (int i = 0; i < chemistry_element_count; i++) {
+ sp->feedback_data.to_distribute.metal_mass[i] += metal_mass_released[i] * norm_factor;
+ }
+ for (int i = 0; i < chemistry_element_count; i++) {
+ sp->feedback_data.to_distribute.mass_from_SNII += sp->feedback_data.to_distribute.metal_mass[i];
+ }
+ sp->feedback_data.to_distribute.total_metal_mass +=
+ metal_mass_released_total * norm_factor;
+ sp->feedback_data.to_distribute.metal_mass_from_SNII +=
+ metal_mass_released_total * norm_factor;
+ } else {
+ error("wrong normalization!!!! mass_released = %e\n", mass_released);
+ }
+}
+
+/**
+ * @brief compute enrichment and feedback due to AGB. To do this, integrate the
+ * IMF weighted by the yields read from tables for each of the quantities of
+ * interest.
+ *
+ * @param log10_min_mass log10 mass at the end of step
+ * @param log10_max_mass log10 mass at the beginning of step
+ * @param stellar_yields array to store calculated yields for passing to
+ * integrate_imf
+ * @param feedback_props star properties data structure
+ * @param sp spart we are computing feedback from
+ */
+inline static void evolve_AGB(float log10_min_mass, float log10_max_mass,
+ float* stellar_yields,
+ const struct feedback_props* restrict feedback_props,
+ struct spart* restrict sp) {
+ int low_imf_mass_bin_index, high_imf_mass_bin_index, mass_bin_index;
+
+ /* If mass at end of step is greater than tabulated lower bound for IMF, limit
+ * it.*/
+ if (log10_max_mass > feedback_props->log10_SNII_min_mass_msun)
+ log10_max_mass = feedback_props->log10_SNII_min_mass_msun;
+
+ /* Don't do anything if the stellar mass hasn't decreased by the end of the
+ * step */
+ if (log10_min_mass >= log10_max_mass) return;
+
+ /* determine which IMF mass bins contribute to the integral */
+ determine_imf_bins(log10_min_mass, log10_max_mass, &low_imf_mass_bin_index,
+ &high_imf_mass_bin_index, feedback_props);
+
+ /* determine which metallicity bin and offset this star belongs to */
+ int iz_low = 0, iz_high = 0, low_index_3d, high_index_3d, low_index_2d, high_index_2d;
+ float dz = 0.f;
+ determine_bin_yield_AGB(&iz_low, &iz_high, &dz,
+ log10(sp->chemistry_data.metal_mass_fraction_total),
+ feedback_props);
+
+ /* compute metals produced */
+ float metal_mass_released[chemistry_element_count], metal_mass_released_total;
+ for (int elem = 0; elem < chemistry_element_count; elem++) {
+ for (mass_bin_index = low_imf_mass_bin_index;
+ mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) {
+ low_index_3d = feedback_row_major_index_3d(
+ iz_low, elem, mass_bin_index, feedback_props->AGB_n_z,
+ chemistry_element_count, feedback_props->n_imf_mass_bins);
+ high_index_3d = feedback_row_major_index_3d(
+ iz_high, elem, mass_bin_index, feedback_props->AGB_n_z,
+ chemistry_element_count, feedback_props->n_imf_mass_bins);
+ low_index_2d = feedback_row_major_index_2d(
+ iz_low, mass_bin_index, feedback_props->AGB_n_z,
+ feedback_props->n_imf_mass_bins);
+ high_index_2d = feedback_row_major_index_2d(
+ iz_high, mass_bin_index, feedback_props->AGB_n_z,
+ feedback_props->n_imf_mass_bins);
+ stellar_yields[mass_bin_index] =
+ (1 - dz) *
+ (feedback_props->yield_AGB.yield_IMF_resampled[low_index_3d] +
+ sp->chemistry_data.metal_mass_fraction[elem] *
+ feedback_props->yield_AGB
+ .ejecta_IMF_resampled[low_index_2d]) +
+ dz * (feedback_props->yield_AGB.yield_IMF_resampled[high_index_3d] +
+ sp->chemistry_data.metal_mass_fraction[elem] *
+ feedback_props->yield_AGB
+ .ejecta_IMF_resampled[high_index_2d]);
+ }
+
+ metal_mass_released[elem] =
+ integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, feedback_props);
+ }
+
+ /* Compute mass produced */
+ for (mass_bin_index = low_imf_mass_bin_index;
+ mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) {
+ low_index_2d = feedback_row_major_index_2d(iz_low, mass_bin_index,
+ feedback_props->AGB_n_z,
+ feedback_props->n_imf_mass_bins);
+ high_index_2d = feedback_row_major_index_2d(
+ iz_high, mass_bin_index, feedback_props->AGB_n_z,
+ feedback_props->n_imf_mass_bins);
+ stellar_yields[mass_bin_index] =
+ (1 - dz) *
+ (feedback_props->yield_AGB
+ .total_metals_IMF_resampled[low_index_2d] +
+ sp->chemistry_data.metal_mass_fraction_total *
+ feedback_props->yield_AGB.ejecta_IMF_resampled[low_index_2d]) +
+ dz *
+ (feedback_props->yield_AGB
+ .total_metals_IMF_resampled[high_index_2d] +
+ sp->chemistry_data.metal_mass_fraction_total *
+ feedback_props->yield_AGB.ejecta_IMF_resampled[high_index_2d]);
+ }
+
+ metal_mass_released_total =
+ integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, feedback_props);
+
+ /* yield normalization */
+ float mass_ejected, mass_released;
+
+ /* zero all negative values */
+ for (int i = 0; i < chemistry_element_count; i++)
+ metal_mass_released[i] = max(metal_mass_released[i], 0.f);
+
+ metal_mass_released_total = max(metal_mass_released_total, 0.f);
+
+ /* compute the total metal mass ejected from the star */
+ for (mass_bin_index = low_imf_mass_bin_index;
+ mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) {
+ low_index_2d = feedback_row_major_index_2d(iz_low, mass_bin_index,
+ feedback_props->AGB_n_z,
+ feedback_props->n_imf_mass_bins);
+ high_index_2d = feedback_row_major_index_2d(
+ iz_high, mass_bin_index, feedback_props->AGB_n_z,
+ feedback_props->n_imf_mass_bins);
+ stellar_yields[mass_bin_index] =
+ (1 - dz) *
+ feedback_props->yield_AGB.ejecta_IMF_resampled[low_index_2d] +
+ dz * feedback_props->yield_AGB.ejecta_IMF_resampled[high_index_2d];
+ }
+
+ mass_ejected =
+ integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, feedback_props);
+
+ /* compute the total mass released */
+ mass_released = metal_mass_released_total +
+ metal_mass_released[chemistry_element_H] +
+ metal_mass_released[chemistry_element_He];
+
+ /* normalize the yields */
+ if (mass_released > 0) {
+ /* Set normalisation factor. Note additional multiplication by the stellar
+ * initial mass as tables are per initial mass */
+ const float norm_factor = sp->mass_init * mass_ejected / mass_released;
+
+ for (int i = 0; i < chemistry_element_count; i++) {
+ sp->feedback_data.to_distribute.metal_mass[i] += metal_mass_released[i] * norm_factor;
+ sp->feedback_data.to_distribute.mass_from_AGB += metal_mass_released[i] * norm_factor;
+ }
+ sp->feedback_data.to_distribute.total_metal_mass +=
+ metal_mass_released_total * norm_factor;
+ sp->feedback_data.to_distribute.metal_mass_from_AGB +=
+ metal_mass_released_total * norm_factor;
+ } else {
+ error("wrong normalization!!!! mass_released = %e\n", mass_released);
+ }
+}
+
+/**
+ * @brief calculates stellar mass in spart that died over the timestep, calls
+ * functions to calculate feedback due to SNIa, SNII and AGB
+ *
+ * @param star_properties feedback_props_props data structure
+ * @param sp spart that we're evolving
+ * @param us unit_system data structure
+ * @param age age of spart at beginning of step
+ * @param dt length of current timestep
+ */
+void compute_stellar_evolution(
+ const struct feedback_props* feedback_props,
+ struct spart* sp, const struct unit_system* us, float age,
+ double dt) {
+
+ /* Allocate temporary array for calculating imf weights */
+ float* stellar_yields;
+ stellar_yields =
+ malloc(feedback_props->n_imf_mass_bins * sizeof(float));
+
+ /* Convert dt and stellar age from internal units to Gyr. */
+ const double Gyr_in_cgs = 3.154e16;
+ double dt_Gyr =
+ dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME) / Gyr_in_cgs;
+ double star_age_Gyr =
+ age * units_cgs_conversion_factor(us, UNIT_CONV_TIME) / Gyr_in_cgs;
+
+ /* calculate mass of stars that has died from the star's birth up to the
+ * beginning and end of timestep */
+ double log10_max_dying_mass_msun = log10(dying_mass_msun(
+ star_age_Gyr, sp->chemistry_data.metal_mass_fraction_total,
+ feedback_props));
+ double log10_min_dying_mass_msun = log10(dying_mass_msun(
+ star_age_Gyr + dt_Gyr, sp->chemistry_data.metal_mass_fraction_total,
+ feedback_props));
+
+ /* Sanity check. Worth investigating if necessary as functions for evaluating
+ * mass of stars dying might be strictly decreasing. */
+ if (log10_min_dying_mass_msun > log10_max_dying_mass_msun)
+ error("min dying mass is greater than max dying mass");
+
+ /* Integration interval is zero - this can happen if minimum and maximum
+ * dying masses are above imf_max_mass_msun. Return without doing any
+ * feedback. */
+ if (log10_min_dying_mass_msun == log10_max_dying_mass_msun) return;
+
+ /* Evolve SNIa, SNII, AGB */
+ evolve_SNIa(log10_min_dying_mass_msun, log10_max_dying_mass_msun,
+ feedback_props, sp, star_age_Gyr, dt_Gyr);
+ evolve_SNII(log10_min_dying_mass_msun, log10_max_dying_mass_msun,
+ stellar_yields, feedback_props, sp);
+ evolve_AGB(log10_min_dying_mass_msun, log10_max_dying_mass_msun,
+ stellar_yields, feedback_props, sp);
+
+ /* Compute the mass to distribute */
+ sp->feedback_data.to_distribute.mass = sp->feedback_data.to_distribute.total_metal_mass +
+ sp->feedback_data.to_distribute.metal_mass[chemistry_element_H] +
+ sp->feedback_data.to_distribute.metal_mass[chemistry_element_He];
+
+ /* Clean up */
+ free(stellar_yields);
+}
+
+/**
+ * @brief Compute number of SN that should go off given the age of the spart
+ *
+ * @param sp spart we're evolving
+ * @param stars_properties stars_props data structure
+ * @param age age of star at the beginning of the step
+ * @param dt length of step
+ */
+float compute_SNe(struct spart* sp,
+ const struct feedback_props* feedback_props,
+ float age, double dt) {
+
+ if (age <= feedback_props->SNII_wind_delay &&
+ age + dt > feedback_props->SNII_wind_delay) {
+
+ return feedback_props->num_SNII_per_msun * sp->mass_init /
+ feedback_props->const_solar_mass;
+ } else {
+ return 0;
+ }
+}
+
+
+/**
+ * @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 feedback_props* feedback_props) {
+
+ /* Set number of elements found in yield tables */
+ feedback_props->SNIa_n_elements = 42;
+ feedback_props->SNII_n_mass = 11;
+ feedback_props->SNII_n_elements = 11;
+ feedback_props->SNII_n_z = 5;
+ feedback_props->AGB_n_mass = 23;
+ feedback_props->AGB_n_elements = 11;
+ feedback_props->AGB_n_z = 3;
+ feedback_props->lifetimes.n_mass = 30;
+ feedback_props->lifetimes.n_z = 6;
+ feedback_props->element_name_length = 15;
+
+ /* Set bounds for imf */
+ feedback_props->log10_SNII_min_mass_msun = 0.77815125f; // log10(6).
+ feedback_props->log10_SNII_max_mass_msun = 2.f; // log10(100).
+ feedback_props->log10_SNIa_max_mass_msun = 0.90308999f; // log10(8).
+
+ /* Yield table filepath */
+ parser_get_param_string(params, "EAGLEFeedback:filename",
+ feedback_props->yield_table_path);
+ parser_get_param_string(params, "EAGLEFeedback:imf_model",
+ feedback_props->IMF_Model);
+
+ /* Initialise IMF */
+ init_imf(feedback_props);
+
+ /* Allocate yield tables */
+ allocate_yield_tables(feedback_props);
+
+ /* Set factors for each element adjusting SNII yield */
+ feedback_props->typeII_factor[0] = 1.f;
+ feedback_props->typeII_factor[1] = 1.f;
+ feedback_props->typeII_factor[2] = 0.5f;
+ feedback_props->typeII_factor[3] = 1.f;
+ feedback_props->typeII_factor[4] = 1.f;
+ feedback_props->typeII_factor[5] = 1.f;
+ feedback_props->typeII_factor[6] = 2.f;
+ feedback_props->typeII_factor[7] = 1.f;
+ feedback_props->typeII_factor[8] = 0.5f;
+
+ /* Read the tables */
+ read_yield_tables(feedback_props);
+
+ /* Set yield_mass_bins array */
+ const float imf_log10_mass_bin_size =
+ (feedback_props->log10_imf_max_mass_msun -
+ feedback_props->log10_imf_min_mass_msun) /
+ (feedback_props->n_imf_mass_bins - 1);
+ for (int i = 0; i < feedback_props->n_imf_mass_bins; i++)
+ feedback_props->yield_mass_bins[i] =
+ imf_log10_mass_bin_size * i + feedback_props->log10_imf_min_mass_msun;
+
+ /* Resample yields from mass bins used in tables to mass bins used in IMF */
+ compute_yields(feedback_props);
+
+ /* Resample ejecta contribution to enrichment from mass bins used in tables to
+ * mass bins used in IMF */
+ compute_ejecta(feedback_props);
+
+ /* Calculate number of type II SN per solar mass
+ * Note: since we are integrating the IMF without weighting it by the yields
+ * pass NULL pointer for stellar_yields array */
+ feedback_props->num_SNII_per_msun =
+ integrate_imf(feedback_props->log10_SNII_min_mass_msun,
+ feedback_props->log10_SNII_max_mass_msun, 0,
+ /*(stellar_yields=)*/ NULL, feedback_props);
+
+ message("initialized stellar feedback");
+}
+
+/**
+ * @brief Initialize the global properties of the feedback scheme.
+ *
+ * By default, takes the values provided by the hydro.
+ *
+ * @param sp The #feedback_properties.
+ * @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.
+ */
+void feedbakc_props_init(struct feedback_props *fp,
+ const struct phys_const *phys_const,
+ const struct unit_system *us,
+ struct swift_params *params,
+ const struct hydro_props *p,
+ const struct cosmology *cosmo) {
+
+ /* Read SNIa timscale */
+ fp->SNIa_timescale_Gyr =
+ parser_get_param_float(params, "EAGLEFeedback:SNIa_timescale_Gyr");
+
+ /* Read the efficiency of producing SNIa */
+ fp->SNIa_efficiency =
+ parser_get_param_float(params, "EAGLEFeedback:SNIa_efficiency");
+
+ /* Are we doing continuous heating? */
+ fp->continuous_heating =
+ parser_get_param_int(params, "EAGLEFeedback:continuous_heating_switch");
+
+ /* Set the delay time before SNII occur */
+ const float Gyr_in_cgs = 1e9 * 365 * 24 * 3600;
+ fp->SNII_wind_delay =
+ parser_get_param_float(params, "EAGLEFeedback:SNII_wind_delay_Gyr") *
+ Gyr_in_cgs / units_cgs_conversion_factor(us, UNIT_CONV_TIME);
+
+ /* Read the temperature change to use in stochastic heating */
+ fp->SNe_deltaT_desired =
+ parser_get_param_float(params, "EAGLEFeedback:SNe_heating_temperature_K");
+ fp->SNe_deltaT_desired /=
+ units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+
+ /* Set ejecta thermal energy */
+ const float ejecta_velocity =
+ 1.0e6 / units_cgs_conversion_factor(
+ us, UNIT_CONV_SPEED); // EAGLE parameter is 10 km/s
+ fp->ejecta_specific_thermal_energy = 0.5 * ejecta_velocity * ejecta_velocity;
+
+ /* Energy released by supernova */
+ fp->total_energy_SNe =
+ 1.0e51 / units_cgs_conversion_factor(us, UNIT_CONV_ENERGY);
+
+ /* Calculate temperature to internal energy conversion factor */
+ fp->temp_to_u_factor =
+ phys_const->const_boltzmann_k /
+ (p->mu_ionised * (hydro_gamma_minus_one)*phys_const->const_proton_mass);
+
+ /* Read birth time to set all stars in ICs to (defaults to -1 to indicate star
+ * present in ICs) */
+ fp->spart_first_init_birth_time = parser_get_opt_param_float(
+ params, "EAGLEFeedback:birth_time_override", -1);
+
+ /* Copy over solar mass */
+ fp->const_solar_mass = phys_const->const_solar_mass;
+}
+
+
diff --git a/src/feedback/EAGLE/feedback.h b/src/feedback/EAGLE/feedback.h
index 6f669956027185a4218f23fd4baa67750dd0a941..d4ea9cb35d100588242343be6738bc87e61c7d9d 100644
--- a/src/feedback/EAGLE/feedback.h
+++ b/src/feedback/EAGLE/feedback.h
@@ -1,3 +1,39 @@
+/*******************************************************************************
+ * 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 <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_FEEDBACK_EAGLE_H
+#define SWIFT_FEEDBACK_EAGLE_H
+
+#include "cosmology.h"
+#include "units.h"
+#include "part.h"
+
+#include "feedback_properties.h"
+
+
+void compute_stellar_evolution(
+ const struct feedback_props* feedback_props,
+ struct spart* sp, const struct unit_system* us, float age,
+ double dt);
+
+float compute_SNe(struct spart* sp,
+ const struct feedback_props* stars_properties,
+ float age, double dt);
/**
* @brief Prepares a s-particle for its feedback interactions
@@ -7,8 +43,8 @@
__attribute__((always_inline)) INLINE static void feedback_init_spart(
struct spart* sp) {
- sp->density_weighted_frac_normalisation_inv = 0.f;
- sp->ngb_mass = 0.f;
+ sp->feedback_data.density_weighted_frac_normalisation_inv = 0.f;
+ sp->feedback_data.ngb_mass = 0.f;
}
/**
@@ -23,8 +59,8 @@ __attribute__((always_inline)) INLINE static void feedback_init_spart(
__attribute__((always_inline)) INLINE static void feedback_first_init_spart(
struct spart* sp, const struct feedback_props* feedback_props) {
- sp->birth_density = -1.f;
- sp->birth_time = feedback_props.spart_first_init_birth_time;
+ //sp->birth_density = -1.f;
+ //sp->birth_time = feedback_props.spart_first_init_birth_time;
feedback_init_spart(sp);
}
@@ -41,566 +77,6 @@ __attribute__((always_inline)) INLINE static void feedback_first_init_spart(
__attribute__((always_inline)) INLINE static void feedback_prepare_spart(
struct spart* sp, const struct feedback_props* feedback_props) {}
-/**
- * @brief determine which metallicity bin star belongs to for AGB, compute bin
- * indices and offsets
- *
- * @param iz_low Pointer to index of metallicity bin to which the star belongs
- * (to be calculated in this function)
- * @param iz_high Pointer to index of metallicity bin to above the star's
- * metallicity (to be calculated in this function)
- * @param dz metallicity bin offset
- * @param log10_metallicity log10 of star metallicity
- * @param star_properties stars_props data structure
- */
-inline static void determine_bin_yield_AGB(
- int* iz_low, int* iz_high, float* dz, float log10_metallicity,
- const struct stars_props* restrict star_properties) {
-
- if (log10_metallicity > log10_min_metallicity) {
- /* Find metallicity bin which contains the star's metallicity */
- int j;
- for (j = 0; j < star_properties->feedback.AGB_n_z - 1 &&
- log10_metallicity >
- star_properties->feedback.yield_AGB.metallicity[j + 1];
- j++)
- ;
- *iz_low = j;
- *iz_high = *iz_low + 1;
-
- /* Compute offset */
- if (log10_metallicity >=
- star_properties->feedback.yield_AGB.metallicity[0] &&
- log10_metallicity <=
- star_properties->feedback.yield_AGB
- .metallicity[star_properties->feedback.AGB_n_z - 1])
- *dz = log10_metallicity -
- star_properties->feedback.yield_AGB.metallicity[*iz_low];
- else
- *dz = 0;
-
- /* Normalize offset */
- float deltaz = star_properties->feedback.yield_AGB.metallicity[*iz_high] -
- star_properties->feedback.yield_AGB.metallicity[*iz_low];
-
- if (deltaz > 0)
- *dz /= deltaz;
- else
- dz = 0;
- } else {
- *iz_low = 0;
- *iz_high = 0;
- *dz = 0;
- }
-}
-
-/**
- * @brief determine which metallicity bin star belongs to for SNII, compute bin
- * indices and offsets
- *
- * @param iz_low Pointer to index of metallicity bin to which the star belongs
- * (to be calculated in this function)
- * @param iz_high Pointer to index of metallicity bin to above the star's
- * metallicity (to be calculated in this function)
- * @param dz metallicity bin offset
- * @param log10_metallicity log10 of star metallicity
- * @param star_properties stars_props data structure
- */
-inline static void determine_bin_yield_SNII(
- int* iz_low, int* iz_high, float* dz, float log10_metallicity,
- const struct stars_props* restrict star_properties) {
-
- if (log10_metallicity > log10_min_metallicity) {
- /* Find metallicity bin which contains the star's metallicity */
- int j;
- for (j = 0; j < star_properties->feedback.SNII_n_z - 1 &&
- log10_metallicity >
- star_properties->feedback.yield_SNII.metallicity[j + 1];
- j++)
- ;
- *iz_low = j;
- *iz_high = *iz_low + 1;
-
- /* Compute offset */
- if (log10_metallicity >=
- star_properties->feedback.yield_SNII.metallicity[0] &&
- log10_metallicity <=
- star_properties->feedback.yield_SNII
- .metallicity[star_properties->feedback.SNII_n_z - 1])
- *dz = log10_metallicity -
- star_properties->feedback.yield_SNII.metallicity[*iz_low];
- else
- *dz = 0;
-
- /* Normalize offset */
- float deltaz = star_properties->feedback.yield_SNII.metallicity[*iz_high] -
- star_properties->feedback.yield_SNII.metallicity[*iz_low];
-
- if (deltaz > 0)
- *dz = *dz / deltaz;
- else
- dz = 0;
- } else {
- *iz_low = 0;
- *iz_high = 0;
- *dz = 0;
- }
-}
-
-/**
- * @brief compute enrichment and feedback due to SNIa. To do this compute the
- * number of SNIa that occur during the timestep, multiply by constants read
- * from tables.
- *
- * @param log10_min_mass log10 mass at the end of step
- * @param log10_max_mass log10 mass at the beginning of step
- * @param stars star properties data structure
- * @param sp spart we are computing feedback from
- * @param star_age_Gyr age of star in Gyr
- * @param dt_Gyr timestep dt in Gyr
- */
-inline static void evolve_SNIa(float log10_min_mass, float log10_max_mass,
- const struct stars_props* restrict stars,
- struct spart* restrict sp, float star_age_Gyr,
- float dt_Gyr) {
-
- /* Check if we're outside the mass range for SNIa */
- if (log10_min_mass >= stars->feedback.log10_SNIa_max_mass_msun) return;
-
- /* If the max mass is outside the mass range update it to be the maximum and
- * use updated values for the star's age and timestep in this function */
- if (log10_max_mass > stars->feedback.log10_SNIa_max_mass_msun) {
- log10_max_mass = stars->feedback.log10_SNIa_max_mass_msun;
- float lifetime_Gyr =
- lifetime_in_Gyr(exp(M_LN10 * stars->feedback.log10_SNIa_max_mass_msun),
- sp->chemistry_data.metal_mass_fraction_total, stars);
- dt_Gyr = star_age_Gyr + dt_Gyr - lifetime_Gyr;
- star_age_Gyr = lifetime_Gyr;
- }
-
- /* compute the number of SNIa */
- /* Efolding (Forster 2006) */
- float num_SNIa_per_msun =
- stars->feedback.SNIa_efficiency *
- (exp(-star_age_Gyr / stars->feedback.SNIa_timescale_Gyr) -
- exp(-(star_age_Gyr + dt_Gyr) / stars->feedback.SNIa_timescale_Gyr)) *
- sp->mass_init;
-
- sp->to_distribute.num_SNIa =
- num_SNIa_per_msun / stars->feedback.const_solar_mass;
-
- /* compute mass fractions of each metal */
- for (int i = 0; i < chemistry_element_count; i++) {
- sp->to_distribute.metal_mass[i] +=
- num_SNIa_per_msun * stars->feedback.yield_SNIa_IMF_resampled[i];
- }
-
- /* Update the metallicity of the material released */
- sp->to_distribute.metal_mass_from_SNIa +=
- num_SNIa_per_msun * stars->feedback.yield_SNIa_total_metals_IMF_resampled;
-
- /* Update the metal mass produced */
- sp->to_distribute.total_metal_mass +=
- num_SNIa_per_msun * stars->feedback.yield_SNIa_total_metals_IMF_resampled;
-
- /* Compute the mass produced by SNIa */
- sp->to_distribute.mass_from_SNIa +=
- num_SNIa_per_msun * stars->feedback.yield_SNIa_total_metals_IMF_resampled;
-
- /* Compute the iron mass produced */
- sp->to_distribute.Fe_mass_from_SNIa +=
- num_SNIa_per_msun *
- stars->feedback.yield_SNIa_IMF_resampled[chemistry_element_Fe];
-}
-
-/**
- * @brief compute enrichment and feedback due to SNII. To do this, integrate the
- * IMF weighted by the yields read from tables for each of the quantities of
- * interest.
- *
- * @param log10_min_mass log10 mass at the end of step
- * @param log10_max_mass log10 mass at the beginning of step
- * @param stellar_yields array to store calculated yields for passing to
- * integrate_imf
- * @param stars star properties data structure
- * @param sp spart we are computing feedback from
- */
-inline static void evolve_SNII(float log10_min_mass, float log10_max_mass,
- float* stellar_yields,
- const struct stars_props* restrict stars,
- struct spart* restrict sp) {
- int low_imf_mass_bin_index, high_imf_mass_bin_index, mass_bin_index;
-
- /* If mass at beginning of step is less than tabulated lower bound for IMF,
- * limit it.*/
- if (log10_min_mass < stars->feedback.log10_SNII_min_mass_msun)
- log10_min_mass = stars->feedback.log10_SNII_min_mass_msun;
-
- /* If mass at end of step is greater than tabulated upper bound for IMF, limit
- * it.*/
- if (log10_max_mass > stars->feedback.log10_SNII_max_mass_msun)
- log10_max_mass = stars->feedback.log10_SNII_max_mass_msun;
-
- /* Don't do anything if the stellar mass hasn't decreased by the end of the
- * step */
- if (log10_min_mass >= log10_max_mass) return;
-
- /* determine which IMF mass bins contribute to the integral */
- determine_imf_bins(log10_min_mass, log10_max_mass, &low_imf_mass_bin_index,
- &high_imf_mass_bin_index, stars);
-
- /* Integrate IMF to determine number of SNII */
- sp->to_distribute.num_SNII =
- integrate_imf(log10_min_mass, log10_max_mass, 0, stellar_yields, stars);
-
- /* determine which metallicity bin and offset this star belongs to */
- int iz_low, iz_high, low_index_3d, high_index_3d, low_index_2d, high_index_2d;
- float dz;
- determine_bin_yield_SNII(&iz_low, &iz_high, &dz,
- log10(sp->chemistry_data.metal_mass_fraction_total),
- stars);
-
- /* compute metals produced */
- float metal_mass_released[chemistry_element_count], metal_mass_released_total;
- for (int elem = 0; elem < chemistry_element_count; elem++) {
- for (mass_bin_index = low_imf_mass_bin_index;
- mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) {
- low_index_3d = feedback_row_major_index_3d(
- iz_low, elem, mass_bin_index, stars->feedback.SNII_n_z,
- chemistry_element_count, stars->feedback.n_imf_mass_bins);
- high_index_3d = feedback_row_major_index_3d(
- iz_high, elem, mass_bin_index, stars->feedback.SNII_n_z,
- chemistry_element_count, stars->feedback.n_imf_mass_bins);
- low_index_2d = feedback_row_major_index_2d(
- iz_low, mass_bin_index, stars->feedback.SNII_n_z,
- stars->feedback.n_imf_mass_bins);
- high_index_2d = feedback_row_major_index_2d(
- iz_high, mass_bin_index, stars->feedback.SNII_n_z,
- stars->feedback.n_imf_mass_bins);
- stellar_yields[mass_bin_index] =
- (1 - dz) *
- (stars->feedback.yield_SNII.yield_IMF_resampled[low_index_3d] +
- sp->chemistry_data.metal_mass_fraction[elem] *
- stars->feedback.yield_SNII
- .ejecta_IMF_resampled[low_index_2d]) +
- dz * (stars->feedback.yield_SNII.yield_IMF_resampled[high_index_3d] +
- sp->chemistry_data.metal_mass_fraction[elem] *
- stars->feedback.yield_SNII
- .ejecta_IMF_resampled[high_index_2d]);
- }
-
- metal_mass_released[elem] =
- integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, stars);
- }
-
- /* Compute mass produced */
- for (mass_bin_index = low_imf_mass_bin_index;
- mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) {
- low_index_2d = feedback_row_major_index_2d(iz_low, mass_bin_index,
- stars->feedback.SNII_n_z,
- stars->feedback.n_imf_mass_bins);
- high_index_2d = feedback_row_major_index_2d(
- iz_high, mass_bin_index, stars->feedback.SNII_n_z,
- stars->feedback.n_imf_mass_bins);
- stellar_yields[mass_bin_index] =
- (1 - dz) * (stars->feedback.yield_SNII
- .total_metals_IMF_resampled[low_index_2d] +
- sp->chemistry_data.metal_mass_fraction_total *
- stars->feedback.yield_SNII
- .ejecta_IMF_resampled[low_index_2d]) +
- dz * (stars->feedback.yield_SNII
- .total_metals_IMF_resampled[high_index_2d] +
- sp->chemistry_data.metal_mass_fraction_total *
- stars->feedback.yield_SNII
- .ejecta_IMF_resampled[high_index_2d]);
- }
-
- metal_mass_released_total =
- integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, stars);
-
- /* yield normalization */
- float mass_ejected, mass_released;
-
- /* zero all negative values */
- for (int i = 0; i < chemistry_element_count; i++)
- metal_mass_released[i] = max(metal_mass_released[i], 0.f);
-
- metal_mass_released_total = max(metal_mass_released_total, 0.f);
-
- /* compute the total metal mass ejected from the star*/
- for (mass_bin_index = low_imf_mass_bin_index;
- mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) {
- low_index_2d = feedback_row_major_index_2d(iz_low, mass_bin_index,
- stars->feedback.SNII_n_z,
- stars->feedback.n_imf_mass_bins);
- high_index_2d = feedback_row_major_index_2d(
- iz_high, mass_bin_index, stars->feedback.SNII_n_z,
- stars->feedback.n_imf_mass_bins);
- stellar_yields[mass_bin_index] =
- (1 - dz) *
- stars->feedback.yield_SNII.ejecta_IMF_resampled[low_index_2d] +
- dz * stars->feedback.yield_SNII.ejecta_IMF_resampled[high_index_2d];
- }
-
- mass_ejected =
- integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, stars);
-
- /* compute the total mass released */
- mass_released = metal_mass_released_total +
- metal_mass_released[chemistry_element_H] +
- metal_mass_released[chemistry_element_He];
-
- /* normalize the yields */
- if (mass_released > 0) {
- /* Set normalisation factor. Note additional multiplication by the star
- * initial mass as tables are per initial mass */
- const float norm_factor = sp->mass_init * mass_ejected / mass_released;
-
- for (int i = 0; i < chemistry_element_count; i++) {
- sp->to_distribute.metal_mass[i] += metal_mass_released[i] * norm_factor;
- }
- for (int i = 0; i < chemistry_element_count; i++) {
- sp->to_distribute.mass_from_SNII += sp->to_distribute.metal_mass[i];
- }
- sp->to_distribute.total_metal_mass +=
- metal_mass_released_total * norm_factor;
- sp->to_distribute.metal_mass_from_SNII +=
- metal_mass_released_total * norm_factor;
- } else {
- error("wrong normalization!!!! mass_released = %e\n", mass_released);
- }
-}
-
-/**
- * @brief compute enrichment and feedback due to AGB. To do this, integrate the
- * IMF weighted by the yields read from tables for each of the quantities of
- * interest.
- *
- * @param log10_min_mass log10 mass at the end of step
- * @param log10_max_mass log10 mass at the beginning of step
- * @param stellar_yields array to store calculated yields for passing to
- * integrate_imf
- * @param stars star properties data structure
- * @param sp spart we are computing feedback from
- */
-inline static void evolve_AGB(float log10_min_mass, float log10_max_mass,
- float* stellar_yields,
- const struct stars_props* restrict stars,
- struct spart* restrict sp) {
- int low_imf_mass_bin_index, high_imf_mass_bin_index, mass_bin_index;
-
- /* If mass at end of step is greater than tabulated lower bound for IMF, limit
- * it.*/
- if (log10_max_mass > stars->feedback.log10_SNII_min_mass_msun)
- log10_max_mass = stars->feedback.log10_SNII_min_mass_msun;
-
- /* Don't do anything if the stellar mass hasn't decreased by the end of the
- * step */
- if (log10_min_mass >= log10_max_mass) return;
-
- /* determine which IMF mass bins contribute to the integral */
- determine_imf_bins(log10_min_mass, log10_max_mass, &low_imf_mass_bin_index,
- &high_imf_mass_bin_index, stars);
-
- /* determine which metallicity bin and offset this star belongs to */
- int iz_low, iz_high, low_index_3d, high_index_3d, low_index_2d, high_index_2d;
- float dz;
- determine_bin_yield_AGB(&iz_low, &iz_high, &dz,
- log10(sp->chemistry_data.metal_mass_fraction_total),
- stars);
-
- /* compute metals produced */
- float metal_mass_released[chemistry_element_count], metal_mass_released_total;
- for (int elem = 0; elem < chemistry_element_count; elem++) {
- for (mass_bin_index = low_imf_mass_bin_index;
- mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) {
- low_index_3d = feedback_row_major_index_3d(
- iz_low, elem, mass_bin_index, stars->feedback.AGB_n_z,
- chemistry_element_count, stars->feedback.n_imf_mass_bins);
- high_index_3d = feedback_row_major_index_3d(
- iz_high, elem, mass_bin_index, stars->feedback.AGB_n_z,
- chemistry_element_count, stars->feedback.n_imf_mass_bins);
- low_index_2d = feedback_row_major_index_2d(
- iz_low, mass_bin_index, stars->feedback.AGB_n_z,
- stars->feedback.n_imf_mass_bins);
- high_index_2d = feedback_row_major_index_2d(
- iz_high, mass_bin_index, stars->feedback.AGB_n_z,
- stars->feedback.n_imf_mass_bins);
- stellar_yields[mass_bin_index] =
- (1 - dz) *
- (stars->feedback.yield_AGB.yield_IMF_resampled[low_index_3d] +
- sp->chemistry_data.metal_mass_fraction[elem] *
- stars->feedback.yield_AGB
- .ejecta_IMF_resampled[low_index_2d]) +
- dz * (stars->feedback.yield_AGB.yield_IMF_resampled[high_index_3d] +
- sp->chemistry_data.metal_mass_fraction[elem] *
- stars->feedback.yield_AGB
- .ejecta_IMF_resampled[high_index_2d]);
- }
-
- metal_mass_released[elem] =
- integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, stars);
- }
-
- /* Compute mass produced */
- for (mass_bin_index = low_imf_mass_bin_index;
- mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) {
- low_index_2d = feedback_row_major_index_2d(iz_low, mass_bin_index,
- stars->feedback.AGB_n_z,
- stars->feedback.n_imf_mass_bins);
- high_index_2d = feedback_row_major_index_2d(
- iz_high, mass_bin_index, stars->feedback.AGB_n_z,
- stars->feedback.n_imf_mass_bins);
- stellar_yields[mass_bin_index] =
- (1 - dz) *
- (stars->feedback.yield_AGB
- .total_metals_IMF_resampled[low_index_2d] +
- sp->chemistry_data.metal_mass_fraction_total *
- stars->feedback.yield_AGB.ejecta_IMF_resampled[low_index_2d]) +
- dz *
- (stars->feedback.yield_AGB
- .total_metals_IMF_resampled[high_index_2d] +
- sp->chemistry_data.metal_mass_fraction_total *
- stars->feedback.yield_AGB.ejecta_IMF_resampled[high_index_2d]);
- }
-
- metal_mass_released_total =
- integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, stars);
-
- /* yield normalization */
- float mass_ejected, mass_released;
-
- /* zero all negative values */
- for (int i = 0; i < chemistry_element_count; i++)
- metal_mass_released[i] = max(metal_mass_released[i], 0.f);
-
- metal_mass_released_total = max(metal_mass_released_total, 0.f);
-
- /* compute the total metal mass ejected from the star */
- for (mass_bin_index = low_imf_mass_bin_index;
- mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) {
- low_index_2d = feedback_row_major_index_2d(iz_low, mass_bin_index,
- stars->feedback.AGB_n_z,
- stars->feedback.n_imf_mass_bins);
- high_index_2d = feedback_row_major_index_2d(
- iz_high, mass_bin_index, stars->feedback.AGB_n_z,
- stars->feedback.n_imf_mass_bins);
- stellar_yields[mass_bin_index] =
- (1 - dz) *
- stars->feedback.yield_AGB.ejecta_IMF_resampled[low_index_2d] +
- dz * stars->feedback.yield_AGB.ejecta_IMF_resampled[high_index_2d];
- }
-
- mass_ejected =
- integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, stars);
-
- /* compute the total mass released */
- mass_released = metal_mass_released_total +
- metal_mass_released[chemistry_element_H] +
- metal_mass_released[chemistry_element_He];
-
- /* normalize the yields */
- if (mass_released > 0) {
- /* Set normalisation factor. Note additional multiplication by the stellar
- * initial mass as tables are per initial mass */
- const float norm_factor = sp->mass_init * mass_ejected / mass_released;
-
- for (int i = 0; i < chemistry_element_count; i++) {
- sp->to_distribute.metal_mass[i] += metal_mass_released[i] * norm_factor;
- sp->to_distribute.mass_from_AGB += metal_mass_released[i] * norm_factor;
- }
- sp->to_distribute.total_metal_mass +=
- metal_mass_released_total * norm_factor;
- sp->to_distribute.metal_mass_from_AGB +=
- metal_mass_released_total * norm_factor;
- } else {
- error("wrong normalization!!!! mass_released = %e\n", mass_released);
- }
-}
-
-/**
- * @brief calculates stellar mass in spart that died over the timestep, calls
- * functions to calculate feedback due to SNIa, SNII and AGB
- *
- * @param star_properties stars_props data structure
- * @param sp spart that we're evolving
- * @param us unit_system data structure
- * @param age age of spart at beginning of step
- * @param dt length of current timestep
- */
-inline static void compute_stellar_evolution(
- const struct stars_props* restrict star_properties,
- struct spart* restrict sp, const struct unit_system* us, float age,
- double dt) {
-
- /* Allocate temporary array for calculating imf weights */
- float* stellar_yields;
- stellar_yields =
- malloc(star_properties->feedback.n_imf_mass_bins * sizeof(float));
-
- /* Convert dt and stellar age from internal units to Gyr. */
- const double Gyr_in_cgs = 3.154e16;
- double dt_Gyr =
- dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME) / Gyr_in_cgs;
- double star_age_Gyr =
- age * units_cgs_conversion_factor(us, UNIT_CONV_TIME) / Gyr_in_cgs;
-
- /* calculate mass of stars that has died from the star's birth up to the
- * beginning and end of timestep */
- double log10_max_dying_mass_msun = log10(dying_mass_msun(
- star_age_Gyr, sp->chemistry_data.metal_mass_fraction_total,
- star_properties));
- double log10_min_dying_mass_msun = log10(dying_mass_msun(
- star_age_Gyr + dt_Gyr, sp->chemistry_data.metal_mass_fraction_total,
- star_properties));
-
- /* Sanity check. Worth investigating if necessary as functions for evaluating
- * mass of stars dying might be strictly decreasing. */
- if (log10_min_dying_mass_msun > log10_max_dying_mass_msun)
- error("min dying mass is greater than max dying mass");
-
- /* Integration interval is zero - this can happen if minimum and maximum
- * dying masses are above imf_max_mass_msun. Return without doing any
- * feedback. */
- if (log10_min_dying_mass_msun == log10_max_dying_mass_msun) return;
-
- /* Evolve SNIa, SNII, AGB */
- evolve_SNIa(log10_min_dying_mass_msun, log10_max_dying_mass_msun,
- star_properties, sp, star_age_Gyr, dt_Gyr);
- evolve_SNII(log10_min_dying_mass_msun, log10_max_dying_mass_msun,
- stellar_yields, star_properties, sp);
- evolve_AGB(log10_min_dying_mass_msun, log10_max_dying_mass_msun,
- stellar_yields, star_properties, sp);
-
- /* Compute the mass to distribute */
- sp->to_distribute.mass = sp->to_distribute.total_metal_mass +
- sp->to_distribute.metal_mass[chemistry_element_H] +
- sp->to_distribute.metal_mass[chemistry_element_He];
-
- /* Clean up */
- free(stellar_yields);
-}
-
-/**
- * @brief Compute number of SN that should go off given the age of the spart
- *
- * @param sp spart we're evolving
- * @param stars_properties stars_props data structure
- * @param age age of star at the beginning of the step
- * @param dt length of step
- */
-inline static float compute_SNe(struct spart* sp,
- const struct stars_props* stars_properties,
- float age, double dt) {
- if (age <= stars_properties->feedback.SNII_wind_delay &&
- age + dt > stars_properties->feedback.SNII_wind_delay) {
- return stars_properties->feedback.num_SNII_per_msun * sp->mass_init /
- stars_properties->feedback.const_solar_mass;
- } else {
- return 0;
- }
-}
-
/**
* @brief Evolve the stellar properties of a #spart.
*
@@ -612,129 +88,52 @@ inline static float compute_SNe(struct spart* sp,
* @param stars_properties The #stars_props
*/
__attribute__((always_inline)) INLINE static void stars_evolve_spart(
- struct spart* restrict sp, const struct stars_props* stars_properties,
+ struct spart* restrict sp, const struct feedback_props* feedback_props,
const struct cosmology* cosmo, const struct unit_system* us,
double star_age, double dt) {
/* Zero the number of SN and amount of mass that is distributed */
- sp->to_distribute.num_SNIa = 0;
- sp->to_distribute.num_SNII = 0;
- sp->to_distribute.mass = 0;
+ sp->feedback_data.to_distribute.num_SNIa = 0;
+ sp->feedback_data.to_distribute.num_SNII = 0;
+ sp->feedback_data.to_distribute.mass = 0;
/* Zero the enrichment quantities */
for (int i = 0; i < chemistry_element_count; i++)
- sp->to_distribute.metal_mass[i] = 0;
- sp->to_distribute.total_metal_mass = 0;
- sp->to_distribute.mass_from_AGB = 0;
- sp->to_distribute.metal_mass_from_AGB = 0;
- sp->to_distribute.mass_from_SNII = 0;
- sp->to_distribute.metal_mass_from_SNII = 0;
- sp->to_distribute.mass_from_SNIa = 0;
- sp->to_distribute.metal_mass_from_SNIa = 0;
- sp->to_distribute.Fe_mass_from_SNIa = 0;
+ sp->feedback_data.to_distribute.metal_mass[i] = 0;
+ sp->feedback_data.to_distribute.total_metal_mass = 0;
+ sp->feedback_data.to_distribute.mass_from_AGB = 0;
+ sp->feedback_data.to_distribute.metal_mass_from_AGB = 0;
+ sp->feedback_data.to_distribute.mass_from_SNII = 0;
+ sp->feedback_data.to_distribute.metal_mass_from_SNII = 0;
+ sp->feedback_data.to_distribute.mass_from_SNIa = 0;
+ sp->feedback_data.to_distribute.metal_mass_from_SNIa = 0;
+ sp->feedback_data.to_distribute.Fe_mass_from_SNIa = 0;
/* Compute amount of enrichment and feedback that needs to be done in this
* step */
- compute_stellar_evolution(stars_properties, sp, us, star_age, dt);
+ compute_stellar_evolution(feedback_props, sp, us, star_age, dt);
/* Decrease star mass by amount of mass distributed to gas neighbours */
- sp->mass -= sp->to_distribute.mass;
+ sp->mass -= sp->feedback_data.to_distribute.mass;
/* Compute the number of type II SNe that went off */
- sp->to_distribute.num_SNe = compute_SNe(sp, stars_properties, star_age, dt);
+ sp->feedback_data.to_distribute.num_SNe = compute_SNe(sp, feedback_props, star_age, dt);
/* Compute energy change due to thermal and kinetic energy of ejecta */
- sp->to_distribute.d_energy =
- sp->to_distribute.mass *
- (stars_properties->feedback.ejecta_specific_thermal_energy +
+ sp->feedback_data.to_distribute.d_energy =
+ sp->feedback_data.to_distribute.mass *
+ (feedback_props->ejecta_specific_thermal_energy +
0.5 * (sp->v[0] * sp->v[0] + sp->v[1] * sp->v[1] + sp->v[2] * sp->v[2]) *
cosmo->a2_inv);
/* Compute probability of heating neighbouring particles */
- if (dt > 0 && sp->ngb_mass > 0)
- sp->to_distribute.heating_probability =
- stars_properties->feedback.total_energy_SNe *
- sp->to_distribute.num_SNe /
- (stars_properties->feedback.temp_to_u_factor *
- stars_properties->feedback.SNe_deltaT_desired * sp->ngb_mass);
+ if (dt > 0 && sp->feedback_data.ngb_mass > 0)
+ sp->feedback_data.to_distribute.heating_probability =
+ feedback_props->total_energy_SNe *
+ sp->feedback_data.to_distribute.num_SNe /
+ (feedback_props->temp_to_u_factor *
+ feedback_props->SNe_deltaT_desired * sp->feedback_data.ngb_mass);
}
-/**
- * @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) {
-
- /* Set number of elements found in yield tables */
- stars->feedback.SNIa_n_elements = 42;
- stars->feedback.SNII_n_mass = 11;
- stars->feedback.SNII_n_elements = 11;
- stars->feedback.SNII_n_z = 5;
- stars->feedback.AGB_n_mass = 23;
- stars->feedback.AGB_n_elements = 11;
- stars->feedback.AGB_n_z = 3;
- stars->feedback.lifetimes.n_mass = 30;
- stars->feedback.lifetimes.n_z = 6;
- stars->feedback.element_name_length = 15;
-
- /* Set bounds for imf */
- stars->feedback.log10_SNII_min_mass_msun = 0.77815125f; // log10(6).
- stars->feedback.log10_SNII_max_mass_msun = 2.f; // log10(100).
- stars->feedback.log10_SNIa_max_mass_msun = 0.90308999f; // log10(8).
-
- /* Yield table filepath */
- parser_get_param_string(params, "EAGLEFeedback:filename",
- stars->feedback.yield_table_path);
- parser_get_param_string(params, "EAGLEFeedback:imf_model",
- stars->feedback.IMF_Model);
-
- /* Initialise IMF */
- init_imf(stars);
-
- /* Allocate yield tables */
- allocate_yield_tables(stars);
- /* Set factors for each element adjusting SNII yield */
- stars->feedback.typeII_factor[0] = 1.f;
- stars->feedback.typeII_factor[1] = 1.f;
- stars->feedback.typeII_factor[2] = 0.5f;
- stars->feedback.typeII_factor[3] = 1.f;
- stars->feedback.typeII_factor[4] = 1.f;
- stars->feedback.typeII_factor[5] = 1.f;
- stars->feedback.typeII_factor[6] = 2.f;
- stars->feedback.typeII_factor[7] = 1.f;
- stars->feedback.typeII_factor[8] = 0.5f;
-
- /* Read the tables */
- read_yield_tables(stars);
-
- /* Set yield_mass_bins array */
- const float imf_log10_mass_bin_size =
- (stars->feedback.log10_imf_max_mass_msun -
- stars->feedback.log10_imf_min_mass_msun) /
- (stars->feedback.n_imf_mass_bins - 1);
- for (int i = 0; i < stars->feedback.n_imf_mass_bins; i++)
- stars->feedback.yield_mass_bins[i] =
- imf_log10_mass_bin_size * i + stars->feedback.log10_imf_min_mass_msun;
-
- /* Resample yields from mass bins used in tables to mass bins used in IMF */
- compute_yields(stars);
-
- /* Resample ejecta contribution to enrichment from mass bins used in tables to
- * mass bins used in IMF */
- compute_ejecta(stars);
-
- /* Calculate number of type II SN per solar mass
- * Note: since we are integrating the IMF without weighting it by the yields
- * pass NULL pointer for stellar_yields array */
- stars->feedback.num_SNII_per_msun =
- integrate_imf(stars->feedback.log10_SNII_min_mass_msun,
- stars->feedback.log10_SNII_max_mass_msun, 0,
- /*(stellar_yields=)*/ NULL, stars);
-
- message("initialized stellar feedback");
-}
+#endif /* SWIFT_FEEDBACK_EAGLE_H */
diff --git a/src/feedback/EAGLE/feedback_iact.h b/src/feedback/EAGLE/feedback_iact.h
index 2122f418d4b95cb843a35b7ed8aeeb7c2c664350..8550ff93545220ee56cee1680a467cdbf20f8ea8 100644
--- a/src/feedback/EAGLE/feedback_iact.h
+++ b/src/feedback/EAGLE/feedback_iact.h
@@ -16,7 +16,7 @@
* @param ti_current Current integer time value
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_stars_density(
+runner_iact_nonsym_feedback_density(
float r2, const float *dx, float hi, float hj, struct spart *restrict si,
const struct part *restrict pj, const struct cosmology *restrict cosmo,
const struct stars_props *restrict stars_properties,
@@ -42,14 +42,14 @@ runner_iact_nonsym_stars_density(
kernel_deval(uj, &wj, &wj_dx);
/* Add mass of pj to neighbour mass of si */
- si->ngb_mass += hydro_get_mass(pj);
+ si->feedback_data.ngb_mass += hydro_get_mass(pj);
/* 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) si->density_weighted_frac_normalisation_inv += wj / rho;
+ if (rho != 0) si->feedback_data.density_weighted_frac_normalisation_inv += wj / rho;
/* Compute contribution to the density */
si->rho_gas += mj * wi;
@@ -72,7 +72,7 @@ runner_iact_nonsym_stars_density(
* generator
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_stars_feedback(
+runner_iact_nonsym_feedback_apply(
float r2, const float *dx, float hi, float hj,
const struct spart *restrict si, struct part *restrict pj,
const struct cosmology *restrict cosmo,
diff --git a/src/feedback/EAGLE/feedback_properties.h b/src/feedback/EAGLE/feedback_properties.h
index 0ccb4c0b33e34c51d633c73d1d5806e821a75f16..0389edee86ec30994a3c81265f70a49a488222d0 100644
--- a/src/feedback/EAGLE/feedback_properties.h
+++ b/src/feedback/EAGLE/feedback_properties.h
@@ -46,7 +46,7 @@ struct lifetime_table {
double **dyingtime;
};
-struct feedback_properties {
+struct feedback_props {
/* Flag to switch between continuous and stochastic heating */
int continuous_heating;
@@ -150,68 +150,3 @@ struct feedback_properties {
float spart_first_init_birth_time;
};
-/**
- * @brief Initialize the global properties of the feedback scheme.
- *
- * By default, takes the values provided by the hydro.
- *
- * @param sp The #feedback_properties.
- * @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.
- */
-INLINE static void feedbakc_props_init(struct feedback_props *fp,
- const struct phys_const *phys_const,
- const struct unit_system *us,
- struct swift_params *params,
- const struct hydro_props *p,
- const struct cosmology *cosmo) {
-
- /* Read SNIa timscale */
- fp->SNIa_timescale_Gyr =
- parser_get_param_float(params, "EAGLEFeedback:SNIa_timescale_Gyr");
-
- /* Read the efficiency of producing SNIa */
- fp->SNIa_efficiency =
- parser_get_param_float(params, "EAGLEFeedback:SNIa_efficiency");
-
- /* Are we doing continuous heating? */
- fp->continuous_heating =
- parser_get_param_int(params, "EAGLEFeedback:continuous_heating_switch");
-
- /* Set the delay time before SNII occur */
- const float Gyr_in_cgs = 1e9 * 365 * 24 * 3600;
- fp->SNII_wind_delay =
- parser_get_param_float(params, "EAGLEFeedback:SNII_wind_delay_Gyr") *
- Gyr_in_cgs / units_cgs_conversion_factor(us, UNIT_CONV_TIME);
-
- /* Read the temperature change to use in stochastic heating */
- fp->SNe_deltaT_desired =
- parser_get_param_float(params, "EAGLEFeedback:SNe_heating_temperature_K");
- fp->SNe_deltaT_desired /=
- units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
-
- /* Set ejecta thermal energy */
- const float ejecta_velocity =
- 1.0e6 / units_cgs_conversion_factor(
- us, UNIT_CONV_SPEED); // EAGLE parameter is 10 km/s
- fp->ejecta_specific_thermal_energy = 0.5 * ejecta_velocity * ejecta_velocity;
-
- /* Energy released by supernova */
- fp->total_energy_SNe =
- 1.0e51 / units_cgs_conversion_factor(us, UNIT_CONV_ENERGY);
-
- /* Calculate temperature to internal energy conversion factor */
- fp->temp_to_u_factor =
- phys_const->const_boltzmann_k /
- (p->mu_ionised * (hydro_gamma_minus_one)*phys_const->const_proton_mass);
-
- /* Read birth time to set all stars in ICs to (defaults to -1 to indicate star
- * present in ICs) */
- fp->spart_first_init_birth_time = parser_get_opt_param_float(
- params, "EAGLEFeedback:birth_time_override", -1);
-
- /* Copy over solar mass */
- fp->const_solar_mass = phys_const->const_solar_mass;
-}
diff --git a/src/feedback/EAGLE/feedback_struct.h b/src/feedback/EAGLE/feedback_struct.h
index 1e0c8e055fb901c7dd78382aca94bbf3a70c6f1d..093bbfc424f507ab4cb9d670954f7cfa37a816a4 100644
--- a/src/feedback/EAGLE/feedback_struct.h
+++ b/src/feedback/EAGLE/feedback_struct.h
@@ -1,7 +1,26 @@
-
-
-
-struct feedback_struct_data {
+/*******************************************************************************
+ * 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 <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_FEEDBACK_STRUCT_EAGLE_H
+#define SWIFT_FEEDBACK_STRUCT_EAGLE_H
+
+
+struct feedback_part_data {
struct {
@@ -58,4 +77,8 @@ struct feedback_struct_data {
/* total mass (unweighted) of neighbouring gas particles */
float ngb_mass;
-}
+};
+
+
+#endif /* SWIFT_FEEDBACK_STRUCT_EAGLE_H */
+
diff --git a/src/feedback/EAGLE/imf.h b/src/feedback/EAGLE/imf.h
index f84cec786c8c196f33688d582ee705ccae5a2ffb..9914adacae2495337ca1406275a9e906fcd98289 100644
--- a/src/feedback/EAGLE/imf.h
+++ b/src/feedback/EAGLE/imf.h
@@ -20,6 +20,9 @@
#define SWIFT_EAGLE_STARS_IMF_H
#include "interpolate.h"
+#include "minmax.h"
+
+#include <string.h>
/**
* @brief the different weightings allowed for the IMF integration
@@ -42,15 +45,15 @@ enum eagle_imf_integration_type {
*/
inline static void determine_imf_bins(
double log10_min_mass, double log10_max_mass, int *i_min, int *i_max,
- const struct stars_props *star_properties) {
+ const struct feedback_props *feedback_props) {
#ifdef SWIFT_DEBUG_CHECKS
if (log10_min_mass > log10_max_mass)
error("Lower bound higher than larger bound.");
#endif
- const int N_bins = star_properties->feedback.n_imf_mass_bins;
- const float *imf_bins_log10 = star_properties->feedback.imf_mass_bin_log10;
+ const int N_bins = feedback_props->n_imf_mass_bins;
+ const float *imf_bins_log10 = feedback_props->imf_mass_bin_log10;
/* Check whether lower mass is within the IMF mass bin range */
log10_min_mass = max(log10_min_mass, imf_bins_log10[0]);
@@ -87,14 +90,14 @@ inline static float integrate_imf(const float log10_min_mass,
const float log10_max_mass,
const enum eagle_imf_integration_type mode,
const float *stellar_yields,
- const struct stars_props *star_properties) {
+ const struct feedback_props *feedback_props) {
/* Pull out some common terms */
- const int N_bins = star_properties->feedback.n_imf_mass_bins;
- const float *imf = star_properties->feedback.imf;
- const float *imf_mass_bin = star_properties->feedback.imf_mass_bin;
+ const int N_bins = feedback_props->n_imf_mass_bins;
+ const float *imf = feedback_props->imf;
+ const float *imf_mass_bin = feedback_props->imf_mass_bin;
const float *imf_mass_bin_log10 =
- star_properties->feedback.imf_mass_bin_log10;
+ feedback_props->imf_mass_bin_log10;
/* IMF mass bin spacing in log10 space. Assumes uniform spacing. */
const float imf_log10_mass_bin_size =
@@ -103,7 +106,7 @@ inline static float integrate_imf(const float log10_min_mass,
/* Determine bins to integrate over based on integration bounds */
int i_min, i_max;
determine_imf_bins(log10_min_mass, log10_max_mass, &i_min, &i_max,
- star_properties);
+ feedback_props);
/* Array for the integrand */
float integrand[N_bins];
@@ -189,82 +192,82 @@ inline static float integrate_imf(const float log10_min_mass,
* table.
*
* @param star_properties #stars_props data structure */
-inline static void init_imf(struct stars_props *restrict star_properties) {
+inline static void init_imf(struct feedback_props *restrict feedback_props) {
/* Define number of imf mass bins */
- star_properties->feedback.n_imf_mass_bins = 200;
+ feedback_props->n_imf_mass_bins = 200;
/* Define max and min imf masses */
- star_properties->feedback.imf_max_mass_msun = 100.f;
- star_properties->feedback.imf_min_mass_msun = 0.1;
- star_properties->feedback.log10_imf_max_mass_msun =
- log10(star_properties->feedback.imf_max_mass_msun);
- star_properties->feedback.log10_imf_min_mass_msun =
- log10(star_properties->feedback.imf_min_mass_msun);
+ feedback_props->imf_max_mass_msun = 100.f;
+ feedback_props->imf_min_mass_msun = 0.1;
+ feedback_props->log10_imf_max_mass_msun =
+ log10(feedback_props->imf_max_mass_msun);
+ feedback_props->log10_imf_min_mass_msun =
+ log10(feedback_props->imf_min_mass_msun);
/* Compute size of mass bins in log10 space */
const float imf_log10_mass_bin_size =
- (star_properties->feedback.log10_imf_max_mass_msun -
- star_properties->feedback.log10_imf_min_mass_msun) /
- (float)(star_properties->feedback.n_imf_mass_bins - 1);
+ (feedback_props->log10_imf_max_mass_msun -
+ feedback_props->log10_imf_min_mass_msun) /
+ (float)(feedback_props->n_imf_mass_bins - 1);
/* Allocate IMF array */
if (swift_memalign(
- "imf-tables", (void **)&star_properties->feedback.imf,
+ "imf-tables", (void **)&feedback_props->imf,
SWIFT_STRUCT_ALIGNMENT,
- star_properties->feedback.n_imf_mass_bins * sizeof(float)) != 0)
+ feedback_props->n_imf_mass_bins * sizeof(float)) != 0)
error("Failed to allocate IMF bins table");
/* Allocate array to store IMF mass bins */
if (swift_memalign(
- "imf-tables", (void **)&star_properties->feedback.imf_mass_bin,
+ "imf-tables", (void **)&feedback_props->imf_mass_bin,
SWIFT_STRUCT_ALIGNMENT,
- star_properties->feedback.n_imf_mass_bins * sizeof(float)) != 0)
+ feedback_props->n_imf_mass_bins * sizeof(float)) != 0)
error("Failed to allocate IMF bins table");
/* Allocate array to store IMF mass bins in log10 space */
if (swift_memalign(
- "imf-tables", (void **)&star_properties->feedback.imf_mass_bin_log10,
+ "imf-tables", (void **)&feedback_props->imf_mass_bin_log10,
SWIFT_STRUCT_ALIGNMENT,
- star_properties->feedback.n_imf_mass_bins * sizeof(float)) != 0)
+ feedback_props->n_imf_mass_bins * sizeof(float)) != 0)
error("Failed to allocate IMF bins table");
/* Set IMF from Chabrier 2003 */
- if (strcmp(star_properties->feedback.IMF_Model, "Chabrier") == 0) {
- for (int i = 0; i < star_properties->feedback.n_imf_mass_bins; i++) {
+ if (strcmp(feedback_props->IMF_Model, "Chabrier") == 0) {
+ for (int i = 0; i < feedback_props->n_imf_mass_bins; i++) {
const float log10_mass_msun =
- star_properties->feedback.log10_imf_min_mass_msun +
+ feedback_props->log10_imf_min_mass_msun +
i * imf_log10_mass_bin_size;
const float mass_msun = exp10f(log10_mass_msun);
if (mass_msun > 1.0)
- star_properties->feedback.imf[i] = 0.237912 * pow(mass_msun, -2.3);
+ feedback_props->imf[i] = 0.237912 * pow(mass_msun, -2.3);
else
- star_properties->feedback.imf[i] =
+ feedback_props->imf[i] =
0.852464 *
exp((log10(mass_msun) - log10(0.079)) *
(log10(mass_msun) - log10(0.079)) / (-2.0 * pow(0.69, 2))) /
mass_msun;
- star_properties->feedback.imf_mass_bin[i] = mass_msun;
- star_properties->feedback.imf_mass_bin_log10[i] = log10_mass_msun;
+ feedback_props->imf_mass_bin[i] = mass_msun;
+ feedback_props->imf_mass_bin_log10[i] = log10_mass_msun;
}
} else {
error("Invalid IMF model %s. Valid models are: 'Chabrier'\n",
- star_properties->feedback.IMF_Model);
+ feedback_props->IMF_Model);
}
/* Normalize the IMF */
const float norm =
- integrate_imf(star_properties->feedback.log10_imf_min_mass_msun,
- star_properties->feedback.log10_imf_max_mass_msun,
+ integrate_imf(feedback_props->log10_imf_min_mass_msun,
+ feedback_props->log10_imf_max_mass_msun,
eagle_imf_integration_mass_weight,
- /*(stellar_yields=)*/ NULL, star_properties);
+ /*(stellar_yields=)*/ NULL, feedback_props);
- for (int i = 0; i < star_properties->feedback.n_imf_mass_bins; i++)
- star_properties->feedback.imf[i] /= norm;
+ for (int i = 0; i < feedback_props->n_imf_mass_bins; i++)
+ feedback_props->imf[i] /= norm;
}
/**
@@ -279,7 +282,7 @@ inline static void init_imf(struct stars_props *restrict star_properties) {
*/
inline static double dying_mass_msun(const float age_Gyr,
const float metallicity,
- const struct stars_props *star_props) {
+ const struct feedback_props *star_props) {
float metallicity_offset, time_offset_low_metallicity = 0,
time_offset_high_metallicity = 0,
@@ -289,104 +292,104 @@ inline static double dying_mass_msun(const float age_Gyr,
time_index_high_metallicity = -1;
if (age_Gyr <= 0) {
- return star_props->feedback.imf_max_mass_msun;
+ return star_props->imf_max_mass_msun;
}
const float log10_age_yr = log10f(age_Gyr * 1.e9);
- if (metallicity <= star_props->feedback.lifetimes.metallicity[0]) {
+ if (metallicity <= star_props->lifetimes.metallicity[0]) {
metallicity_index = 0;
metallicity_offset = 0.0;
} else if (metallicity >=
- star_props->feedback.lifetimes
- .metallicity[star_props->feedback.lifetimes.n_z - 1]) {
- metallicity_index = star_props->feedback.lifetimes.n_z - 2;
+ star_props->lifetimes
+ .metallicity[star_props->lifetimes.n_z - 1]) {
+ metallicity_index = star_props->lifetimes.n_z - 2;
metallicity_offset = 1.0;
} else {
for (metallicity_index = 0;
- metallicity_index < star_props->feedback.lifetimes.n_z - 1;
+ metallicity_index < star_props->lifetimes.n_z - 1;
metallicity_index++)
- if (star_props->feedback.lifetimes.metallicity[metallicity_index + 1] >
+ if (star_props->lifetimes.metallicity[metallicity_index + 1] >
metallicity)
break;
metallicity_offset =
(metallicity -
- star_props->feedback.lifetimes.metallicity[metallicity_index]) /
- (star_props->feedback.lifetimes.metallicity[metallicity_index + 1] -
- star_props->feedback.lifetimes.metallicity[metallicity_index]);
+ star_props->lifetimes.metallicity[metallicity_index]) /
+ (star_props->lifetimes.metallicity[metallicity_index + 1] -
+ star_props->lifetimes.metallicity[metallicity_index]);
}
if (log10_age_yr >=
- star_props->feedback.lifetimes.dyingtime[metallicity_index][0]) {
+ star_props->lifetimes.dyingtime[metallicity_index][0]) {
time_index_low_metallicity = 0;
time_offset_low_metallicity = 0.0;
} else if (log10_age_yr <=
- star_props->feedback.lifetimes
+ star_props->lifetimes
.dyingtime[metallicity_index]
- [star_props->feedback.lifetimes.n_mass - 1]) {
- time_index_low_metallicity = star_props->feedback.lifetimes.n_mass - 2;
+ [star_props->lifetimes.n_mass - 1]) {
+ time_index_low_metallicity = star_props->lifetimes.n_mass - 2;
time_offset_low_metallicity = 1.0;
}
if (log10_age_yr >=
- star_props->feedback.lifetimes.dyingtime[metallicity_index + 1][0]) {
+ star_props->lifetimes.dyingtime[metallicity_index + 1][0]) {
time_index_high_metallicity = 0;
time_offset_high_metallicity = 0.0;
} else if (log10_age_yr <=
- star_props->feedback.lifetimes
+ star_props->lifetimes
.dyingtime[metallicity_index + 1]
- [star_props->feedback.lifetimes.n_mass - 1]) {
- time_index_high_metallicity = star_props->feedback.lifetimes.n_mass - 2;
+ [star_props->lifetimes.n_mass - 1]) {
+ time_index_high_metallicity = star_props->lifetimes.n_mass - 2;
time_offset_high_metallicity = 1.0;
}
- int i = star_props->feedback.lifetimes.n_mass;
+ int i = star_props->lifetimes.n_mass;
while (i >= 0 && (time_index_low_metallicity == -1 ||
time_index_high_metallicity == -1)) {
i--;
- if (star_props->feedback.lifetimes.dyingtime[metallicity_index][i] >=
+ if (star_props->lifetimes.dyingtime[metallicity_index][i] >=
log10_age_yr &&
time_index_low_metallicity == -1) {
time_index_low_metallicity = i;
time_offset_low_metallicity =
(log10_age_yr -
- star_props->feedback.lifetimes
+ star_props->lifetimes
.dyingtime[metallicity_index][time_index_low_metallicity]) /
- (star_props->feedback.lifetimes
+ (star_props->lifetimes
.dyingtime[metallicity_index][time_index_low_metallicity + 1] -
- star_props->feedback.lifetimes
+ star_props->lifetimes
.dyingtime[metallicity_index][time_index_low_metallicity]);
}
- if (star_props->feedback.lifetimes.dyingtime[metallicity_index + 1][i] >=
+ if (star_props->lifetimes.dyingtime[metallicity_index + 1][i] >=
log10_age_yr &&
time_index_high_metallicity == -1) {
time_index_high_metallicity = i;
time_offset_high_metallicity =
(log10_age_yr -
- star_props->feedback.lifetimes
+ star_props->lifetimes
.dyingtime[metallicity_index + 1][time_index_high_metallicity]) /
- (star_props->feedback.lifetimes
+ (star_props->lifetimes
.dyingtime[metallicity_index + 1]
[time_index_high_metallicity + 1] -
- star_props->feedback.lifetimes
+ star_props->lifetimes
.dyingtime[metallicity_index + 1][time_index_high_metallicity]);
}
}
mass_low_metallicity =
- interpolate_1d(star_props->feedback.lifetimes.mass,
+ interpolate_1d(star_props->lifetimes.mass,
time_index_low_metallicity, time_offset_low_metallicity);
mass_high_metallicity =
- interpolate_1d(star_props->feedback.lifetimes.mass,
+ interpolate_1d(star_props->lifetimes.mass,
time_index_high_metallicity, time_offset_high_metallicity);
float mass = (1.0 - metallicity_offset) * mass_low_metallicity +
metallicity_offset * mass_high_metallicity;
/* Check that we haven't killed too many stars */
- if (mass > star_props->feedback.imf_max_mass_msun)
- mass = star_props->feedback.imf_max_mass_msun;
+ if (mass > star_props->imf_max_mass_msun)
+ mass = star_props->imf_max_mass_msun;
return mass;
}
@@ -400,61 +403,61 @@ inline static double dying_mass_msun(const float age_Gyr,
* @param star_properties the #stars_props data struct
*/
inline static float lifetime_in_Gyr(float mass, float metallicity,
- const struct stars_props *restrict
- star_properties) {
+ const struct feedback_props *restrict
+ feedback_props) {
double time_Gyr = 0, mass_offset, metallicity_offset;
int mass_index, metallicity_index;
- if (mass <= star_properties->feedback.lifetimes.mass[0]) {
+ if (mass <= feedback_props->lifetimes.mass[0]) {
mass_index = 0;
mass_offset = 0.0;
} else if (mass >=
- star_properties->feedback.lifetimes
- .mass[star_properties->feedback.lifetimes.n_mass - 1]) {
- mass_index = star_properties->feedback.lifetimes.n_mass - 2;
+ feedback_props->lifetimes
+ .mass[feedback_props->lifetimes.n_mass - 1]) {
+ mass_index = feedback_props->lifetimes.n_mass - 2;
mass_offset = 1.0;
} else {
for (mass_index = 0;
- mass_index < star_properties->feedback.lifetimes.n_mass - 1;
+ mass_index < feedback_props->lifetimes.n_mass - 1;
mass_index++)
- if (star_properties->feedback.lifetimes.mass[mass_index + 1] > mass)
+ if (feedback_props->lifetimes.mass[mass_index + 1] > mass)
break;
mass_offset =
- (mass - star_properties->feedback.lifetimes.mass[mass_index]) /
- (star_properties->feedback.lifetimes.mass[mass_index + 1] -
- star_properties->feedback.lifetimes.mass[mass_index]);
+ (mass - feedback_props->lifetimes.mass[mass_index]) /
+ (feedback_props->lifetimes.mass[mass_index + 1] -
+ feedback_props->lifetimes.mass[mass_index]);
}
- if (metallicity <= star_properties->feedback.lifetimes.metallicity[0]) {
+ if (metallicity <= feedback_props->lifetimes.metallicity[0]) {
metallicity_index = 0;
metallicity_offset = 0.0;
} else if (metallicity >=
- star_properties->feedback.lifetimes
- .metallicity[star_properties->feedback.lifetimes.n_z - 1]) {
- metallicity_index = star_properties->feedback.lifetimes.n_z - 2;
+ feedback_props->lifetimes
+ .metallicity[feedback_props->lifetimes.n_z - 1]) {
+ metallicity_index = feedback_props->lifetimes.n_z - 2;
metallicity_offset = 1.0;
} else {
for (metallicity_index = 0;
- metallicity_index < star_properties->feedback.lifetimes.n_z - 1;
+ metallicity_index < feedback_props->lifetimes.n_z - 1;
metallicity_index++)
- if (star_properties->feedback.lifetimes
+ if (feedback_props->lifetimes
.metallicity[metallicity_index + 1] > metallicity)
break;
metallicity_offset =
(metallicity -
- star_properties->feedback.lifetimes.metallicity[metallicity_index]) /
- (star_properties->feedback.lifetimes
+ feedback_props->lifetimes.metallicity[metallicity_index]) /
+ (feedback_props->lifetimes
.metallicity[metallicity_index + 1] -
- star_properties->feedback.lifetimes.metallicity[metallicity_index]);
+ feedback_props->lifetimes.metallicity[metallicity_index]);
}
/* time in Gyr */
time_Gyr =
- exp(M_LN10 * interpolate_2d(star_properties->feedback.lifetimes.dyingtime,
+ exp(M_LN10 * interpolate_2d(feedback_props->lifetimes.dyingtime,
metallicity_index, mass_index,
metallicity_offset, mass_offset)) /
1.0e9;
diff --git a/src/feedback/EAGLE/interpolate.h b/src/feedback/EAGLE/interpolate.h
index e544d8e2e7cafc2b9e655cd2e2497c720fce14bd..10cf65454c8e3ebed465be885fc721038f7e2379 100644
--- a/src/feedback/EAGLE/interpolate.h
+++ b/src/feedback/EAGLE/interpolate.h
@@ -16,8 +16,40 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_EAGLE_STARS_INTERPOLATE_H
-#define SWIFT_EAGLE_STARS_INTERPOLATE_H
+#ifndef SWIFT_EAGLE_FEEDBACK_INTERPOLATE_H
+#define SWIFT_EAGLE_FEEDBACK_INTERPOLATE_H
+
+#include "error.h"
+
+
+/**
+ * @brief Returns the 1d index of element with 2d indices i,j
+ * from a flattened 2d array in row major order
+ *
+ * ALEXEI: this function also appears in EAGLE cooling. Could this be done
+ * without duplication?
+ *
+ * @param i, j Indices of element of interest
+ * @param nx, ny Sizes of array dimensions
+ */
+__attribute__((always_inline)) INLINE int feedback_row_major_index_2d(int i,
+ int j,
+ int nx,
+ int ny) {
+ return i * ny + j;
+}
+
+/**
+ * @brief Returns the 1d index of element with 3d indices i,j,k
+ * from a flattened 3d array in row major order
+ *
+ * @param i, j, k Indices of element of interest
+ * @param nx, ny, nz Sizes of array dimensions
+ */
+__attribute__((always_inline)) INLINE int feedback_row_major_index_3d(
+ int i, int j, int k, int nx, int ny, int nz) {
+ return i * ny * nz + j * nz + k;
+}
/**
* @brief linear interpolation of 1d table at bin i with offset dx
@@ -94,4 +126,4 @@ inline static double interpolate_1D_non_uniform(double* array_x,
return result;
}
-#endif /* SWIFT_EAGLE_STARS_INTERPOLATE_H */
+#endif /* SWIFT_EAGLE_FEEDBACK_INTERPOLATE_H */
diff --git a/src/feedback/EAGLE/yield_tables.h b/src/feedback/EAGLE/yield_tables.h
index c82868572fd418a0d2e65ee2d71ed89a597d722b..5844949d8155887ecf1f534d9d71726816e07136 100644
--- a/src/feedback/EAGLE/yield_tables.h
+++ b/src/feedback/EAGLE/yield_tables.h
@@ -23,35 +23,6 @@
static const float log10_min_metallicity = -20;
-/**
- * @brief Returns the 1d index of element with 2d indices i,j
- * from a flattened 2d array in row major order
- *
- * ALEXEI: this function also appears in EAGLE cooling. Could this be done
- * without duplication?
- *
- * @param i, j Indices of element of interest
- * @param nx, ny Sizes of array dimensions
- */
-__attribute__((always_inline)) INLINE int feedback_row_major_index_2d(int i,
- int j,
- int nx,
- int ny) {
- return i * ny + j;
-}
-
-/**
- * @brief Returns the 1d index of element with 3d indices i,j,k
- * from a flattened 3d array in row major order
- *
- * @param i, j, k Indices of element of interest
- * @param nx, ny, nz Sizes of array dimensions
- */
-__attribute__((always_inline)) INLINE int feedback_row_major_index_3d(
- int i, int j, int k, int nx, int ny, int nz) {
- return i * ny * nz + j * nz + k;
-}
-
/**
* @brief returns index of element_name within array of element names
* (element_array)
@@ -78,7 +49,8 @@ inline static int get_element_index(const char *element_name,
*
* @param stars the #stars_props data structure
*/
-inline static void read_yield_tables(struct stars_props *restrict stars) {
+inline static void read_yield_tables(struct feedback_props *feedback_props) {
+
#ifdef HAVE_HDF5
int i, j, k, flat_index;
@@ -90,7 +62,7 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
herr_t status;
/* Open SNIa tables for reading */
- sprintf(fname, "%s/SNIa.hdf5", stars->feedback.yield_table_path);
+ sprintf(fname, "%s/SNIa.hdf5", feedback_props->yield_table_path);
file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
if (file_id < 0) error("unable to open file %s\n", fname);
@@ -99,7 +71,7 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
H5Tset_size(datatype, H5T_VARIABLE);
dataset = H5Dopen(file_id, "Species_names", H5P_DEFAULT);
status = H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->feedback.SNIa_element_names);
+ feedback_props->SNIa_element_names);
if (status < 0) error("error reading SNIa element names");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -109,7 +81,7 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
/* read SNIa yields */
dataset = H5Dopen(file_id, "Yield", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->feedback.yields_SNIa);
+ feedback_props->yields_SNIa);
if (status < 0) error("error reading SNIa yields");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -117,7 +89,7 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
/* read SNIa total metals released */
dataset = H5Dopen(file_id, "Total_Metals", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- &stars->feedback.yield_SNIa_total_metals_IMF_resampled);
+ &feedback_props->yield_SNIa_total_metals_IMF_resampled);
if (status < 0) error("error reading SNIa total metal");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -126,7 +98,7 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
if (status < 0) error("error closing SNIa file");
/* Open SNII tables for reading */
- sprintf(fname, "%s/SNII.hdf5", stars->feedback.yield_table_path);
+ sprintf(fname, "%s/SNII.hdf5", feedback_props->yield_table_path);
file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
if (file_id < 0) error("unable to open file %s\n", fname);
@@ -135,7 +107,7 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
H5Tset_size(datatype, H5T_VARIABLE);
dataset = H5Dopen(file_id, "Species_names", H5P_DEFAULT);
status = H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->feedback.SNII_element_names);
+ feedback_props->SNII_element_names);
if (status < 0) error("error reading SNII element names");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -145,7 +117,7 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
/* read array of masses */
dataset = H5Dopen(file_id, "Masses", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->feedback.yield_SNII.mass);
+ feedback_props->yield_SNII.mass);
if (status < 0) error("error reading SNII masses");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -153,17 +125,17 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
/* read array of metallicities */
dataset = H5Dopen(file_id, "Metallicities", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->feedback.yield_SNII.metallicity);
+ feedback_props->yield_SNII.metallicity);
if (status < 0) error("error reading SNII metallicities");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
/* declare temporary arrays to read data from HDF5 files */
- double temp_yield_SNII[stars->feedback.SNII_n_elements]
- [stars->feedback.SNII_n_mass];
- double temp_ejecta_SNII[stars->feedback.SNII_n_mass],
- tempmet1[stars->feedback.SNII_n_mass];
- char *metallicity_yield_table_name_SNII[stars->feedback.SNII_n_z];
+ double temp_yield_SNII[feedback_props->SNII_n_elements]
+ [feedback_props->SNII_n_mass];
+ double temp_ejecta_SNII[feedback_props->SNII_n_mass],
+ tempmet1[feedback_props->SNII_n_mass];
+ char *metallicity_yield_table_name_SNII[feedback_props->SNII_n_z];
/* read metallicity names */
datatype = H5Tcopy(H5T_C_S1);
@@ -178,7 +150,7 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
if (status < 0) error("error closing datatype");
/* read SNII yield tables */
- for (i = 0; i < stars->feedback.SNII_n_z; i++) {
+ for (i = 0; i < feedback_props->SNII_n_z; i++) {
/* read yields to temporary array */
sprintf(setname, "/Yields/%s/Yield", metallicity_yield_table_name_SNII[i]);
dataset = H5Dopen(file_id, setname, H5P_DEFAULT);
@@ -209,17 +181,17 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
if (status < 0) error("error closing dataset");
/* Flatten the temporary tables that were read, store in stars_props */
- for (k = 0; k < stars->feedback.SNII_n_mass; k++) {
- flat_index = feedback_row_major_index_2d(i, k, stars->feedback.SNII_n_z,
- stars->feedback.SNII_n_mass);
- stars->feedback.yield_SNII.ejecta[flat_index] = temp_ejecta_SNII[k];
- stars->feedback.yield_SNII.total_metals[flat_index] = tempmet1[k];
+ for (k = 0; k < feedback_props->SNII_n_mass; k++) {
+ flat_index = feedback_row_major_index_2d(i, k, feedback_props->SNII_n_z,
+ feedback_props->SNII_n_mass);
+ feedback_props->yield_SNII.ejecta[flat_index] = temp_ejecta_SNII[k];
+ feedback_props->yield_SNII.total_metals[flat_index] = tempmet1[k];
- for (j = 0; j < stars->feedback.SNII_n_elements; j++) {
+ for (j = 0; j < feedback_props->SNII_n_elements; j++) {
flat_index = feedback_row_major_index_3d(
- i, j, k, stars->feedback.SNII_n_z, stars->feedback.SNII_n_elements,
- stars->feedback.SNII_n_mass);
- stars->feedback.yield_SNII.yield[flat_index] = temp_yield_SNII[j][k];
+ i, j, k, feedback_props->SNII_n_z, feedback_props->SNII_n_elements,
+ feedback_props->SNII_n_mass);
+ feedback_props->yield_SNII.yield[flat_index] = temp_yield_SNII[j][k];
}
}
}
@@ -228,7 +200,7 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
if (status < 0) error("error closing file");
/* Read AGB tables */
- sprintf(fname, "%s/AGB.hdf5", stars->feedback.yield_table_path);
+ sprintf(fname, "%s/AGB.hdf5", feedback_props->yield_table_path);
file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
if (file_id < 0) error("unable to open file %s\n", fname);
@@ -237,7 +209,7 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
H5Tset_size(datatype, H5T_VARIABLE);
dataset = H5Dopen(file_id, "Species_names", H5P_DEFAULT);
status = H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->feedback.AGB_element_names);
+ feedback_props->AGB_element_names);
if (status < 0) error("error reading AGB element names");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -247,7 +219,7 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
/* read array of masses */
dataset = H5Dopen(file_id, "Masses", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->feedback.yield_AGB.mass);
+ feedback_props->yield_AGB.mass);
if (status < 0) error("error reading AGB masses");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -255,17 +227,17 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
/* read array of metallicities */
dataset = H5Dopen(file_id, "Metallicities", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->feedback.yield_AGB.metallicity);
+ feedback_props->yield_AGB.metallicity);
if (status < 0) error("error reading AGB metallicities");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
/* declare temporary arrays to read data from HDF5 files */
- double temp_yield_AGB[stars->feedback.AGB_n_elements]
- [stars->feedback.AGB_n_mass];
- double temp_ejecta_AGB[stars->feedback.AGB_n_mass],
- tempmet2[stars->feedback.AGB_n_mass];
- char *metallicity_yield_table_name_AGB[stars->feedback.AGB_n_z];
+ double temp_yield_AGB[feedback_props->AGB_n_elements]
+ [feedback_props->AGB_n_mass];
+ double temp_ejecta_AGB[feedback_props->AGB_n_mass],
+ tempmet2[feedback_props->AGB_n_mass];
+ char *metallicity_yield_table_name_AGB[feedback_props->AGB_n_z];
/* read metallicity names */
datatype = H5Tcopy(H5T_C_S1);
@@ -280,7 +252,7 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
if (status < 0) error("error closing datatype");
/* read AGB yield tables */
- for (i = 0; i < stars->feedback.AGB_n_z; i++) {
+ for (i = 0; i < feedback_props->AGB_n_z; i++) {
/* read yields to temporary array */
sprintf(setname, "/Yields/%s/Yield", metallicity_yield_table_name_AGB[i]);
dataset = H5Dopen(file_id, setname, H5P_DEFAULT);
@@ -311,17 +283,17 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
if (status < 0) error("error closing dataset");
/* Flatten the temporary tables that were read, store in stars_props */
- for (k = 0; k < stars->feedback.AGB_n_mass; k++) {
- flat_index = feedback_row_major_index_2d(i, k, stars->feedback.AGB_n_z,
- stars->feedback.AGB_n_mass);
- stars->feedback.yield_AGB.ejecta[flat_index] = temp_ejecta_AGB[k];
- stars->feedback.yield_AGB.total_metals[flat_index] = tempmet2[k];
+ for (k = 0; k < feedback_props->AGB_n_mass; k++) {
+ flat_index = feedback_row_major_index_2d(i, k, feedback_props->AGB_n_z,
+ feedback_props->AGB_n_mass);
+ feedback_props->yield_AGB.ejecta[flat_index] = temp_ejecta_AGB[k];
+ feedback_props->yield_AGB.total_metals[flat_index] = tempmet2[k];
- for (j = 0; j < stars->feedback.AGB_n_elements; j++) {
+ for (j = 0; j < feedback_props->AGB_n_elements; j++) {
flat_index = feedback_row_major_index_3d(
- i, j, k, stars->feedback.AGB_n_z, stars->feedback.AGB_n_elements,
- stars->feedback.AGB_n_mass);
- stars->feedback.yield_AGB.yield[flat_index] = temp_yield_AGB[j][k];
+ i, j, k, feedback_props->AGB_n_z, feedback_props->AGB_n_elements,
+ feedback_props->AGB_n_mass);
+ feedback_props->yield_AGB.yield[flat_index] = temp_yield_AGB[j][k];
}
}
}
@@ -330,14 +302,14 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
if (status < 0) error("error closing file");
/* open lifetimes table */
- sprintf(fname, "%s/Lifetimes.hdf5", stars->feedback.yield_table_path);
+ sprintf(fname, "%s/Lifetimes.hdf5", feedback_props->yield_table_path);
file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
if (file_id < 0) error("unable to open file %s\n", fname);
/* read lifetimes mass bins */
dataset = H5Dopen(file_id, "Masses", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->feedback.lifetimes.mass);
+ feedback_props->lifetimes.mass);
if (status < 0) error("error reading lifetime table masses");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -345,23 +317,23 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
/* read metallicity bins */
dataset = H5Dopen(file_id, "Metallicities", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->feedback.lifetimes.metallicity);
+ feedback_props->lifetimes.metallicity);
if (status < 0) error("error reading lifetimes metallicities");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
/* allocate temporary array to read lifetimes */
- double temp_lifetimes[stars->feedback.lifetimes.n_z]
- [stars->feedback.lifetimes.n_mass];
+ double temp_lifetimes[feedback_props->lifetimes.n_z]
+ [feedback_props->lifetimes.n_mass];
dataset = H5Dopen(file_id, "Lifetimes", H5P_DEFAULT);
H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
temp_lifetimes);
H5Dclose(dataset);
- for (i = 0; i < stars->feedback.lifetimes.n_z; i++) {
- for (j = 0; j < stars->feedback.lifetimes.n_mass; j++) {
- stars->feedback.lifetimes.dyingtime[i][j] = log10(temp_lifetimes[i][j]);
+ for (i = 0; i < feedback_props->lifetimes.n_z; i++) {
+ for (j = 0; j < feedback_props->lifetimes.n_mass; j++) {
+ feedback_props->lifetimes.dyingtime[i][j] = log10(temp_lifetimes[i][j]);
}
}
@@ -375,80 +347,80 @@ inline static void read_yield_tables(struct stars_props *restrict stars) {
*
* @param stars the #stars_props data struct to store the tables in
*/
-inline static void allocate_yield_tables(struct stars_props *restrict stars) {
+inline static void allocate_yield_tables(struct feedback_props *feedback_props) {
/* Allocate array to store SNIa yield tables */
- if (swift_memalign("feedback-tables", (void **)&stars->feedback.yields_SNIa,
+ if (swift_memalign("feedback-tables", (void **)&feedback_props->yields_SNIa,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.SNIa_n_elements * sizeof(double)) != 0) {
+ feedback_props->SNIa_n_elements * sizeof(double)) != 0) {
error("Failed to allocate SNIa yields array");
}
/* Allocate array to store SNIa yield table resampled by IMF mass bins */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_SNIa_IMF_resampled,
+ (void **)&feedback_props->yield_SNIa_IMF_resampled,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.SNIa_n_elements * sizeof(double)) != 0) {
+ feedback_props->SNIa_n_elements * sizeof(double)) != 0) {
error("Failed to allocate SNIa IMF resampled yields array");
}
/* Allocate array for AGB mass bins */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_AGB.mass,
+ (void **)&feedback_props->yield_AGB.mass,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.AGB_n_mass * sizeof(double)) != 0) {
+ feedback_props->AGB_n_mass * sizeof(double)) != 0) {
error("Failed to allocate AGB mass array");
}
/* Allocate array for AGB metallicity bins */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_AGB.metallicity,
+ (void **)&feedback_props->yield_AGB.metallicity,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.AGB_n_z * sizeof(double)) != 0) {
+ feedback_props->AGB_n_z * sizeof(double)) != 0) {
error("Failed to allocate AGB metallicity array");
}
/* Allocate array to store AGB yield tables */
if (swift_memalign(
- "feedback-tables", (void **)&stars->feedback.yield_AGB.yield,
+ "feedback-tables", (void **)&feedback_props->yield_AGB.yield,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.AGB_n_z * stars->feedback.AGB_n_mass *
- stars->feedback.AGB_n_elements * sizeof(double)) != 0) {
+ feedback_props->AGB_n_z * feedback_props->AGB_n_mass *
+ feedback_props->AGB_n_elements * sizeof(double)) != 0) {
error("Failed to allocate AGB yield array");
}
/* Allocate array to store AGB yield table resampled by IMF mass bins */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_AGB.yield_IMF_resampled,
+ (void **)&feedback_props->yield_AGB.yield_IMF_resampled,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.AGB_n_z * stars->feedback.n_imf_mass_bins *
+ feedback_props->AGB_n_z * feedback_props->n_imf_mass_bins *
chemistry_element_count * sizeof(double)) != 0) {
error("Failed to allocate AGB IMF resampled array");
}
/* Allocate array to store AGB ejecta tables */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_AGB.ejecta,
+ (void **)&feedback_props->yield_AGB.ejecta,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.AGB_n_z * stars->feedback.AGB_n_mass *
+ feedback_props->AGB_n_z * feedback_props->AGB_n_mass *
sizeof(double)) != 0) {
error("Failed to allocate AGB ejecta array");
}
/* Allocate array to store AGB ejecta table resampled by IMF mass bins */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_AGB.ejecta_IMF_resampled,
+ (void **)&feedback_props->yield_AGB.ejecta_IMF_resampled,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.AGB_n_z * stars->feedback.n_imf_mass_bins *
+ feedback_props->AGB_n_z * feedback_props->n_imf_mass_bins *
sizeof(double)) != 0) {
error("Failed to allocate AGB ejecta IMF resampled array");
}
/* Allocate array to store table of total metals released by AGB */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_AGB.total_metals,
+ (void **)&feedback_props->yield_AGB.total_metals,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.AGB_n_z * stars->feedback.AGB_n_mass *
+ feedback_props->AGB_n_z * feedback_props->AGB_n_mass *
sizeof(double)) != 0) {
error("Failed to allocate AGB total metals array");
}
@@ -457,72 +429,72 @@ inline static void allocate_yield_tables(struct stars_props *restrict stars) {
* IMF mass bins */
if (swift_memalign(
"feedback-tables",
- (void **)&stars->feedback.yield_AGB.total_metals_IMF_resampled,
+ (void **)&feedback_props->yield_AGB.total_metals_IMF_resampled,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.AGB_n_z * stars->feedback.n_imf_mass_bins *
+ feedback_props->AGB_n_z * feedback_props->n_imf_mass_bins *
sizeof(double)) != 0) {
error("Failed to allocate AGB total metals IMF resampled array");
}
/* Allocate array for SNII mass bins */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_SNII.mass,
+ (void **)&feedback_props->yield_SNII.mass,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.SNII_n_mass * sizeof(double)) != 0) {
+ feedback_props->SNII_n_mass * sizeof(double)) != 0) {
error("Failed to allocate SNII mass array");
}
/* Allocate array for SNII metallicity bins */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_SNII.metallicity,
+ (void **)&feedback_props->yield_SNII.metallicity,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.SNII_n_z * sizeof(double)) != 0) {
+ feedback_props->SNII_n_z * sizeof(double)) != 0) {
error("Failed to allocate SNII metallicity array");
}
/* Allocate array to store SNII yield tables */
if (swift_memalign(
- "feedback-tables", (void **)&stars->feedback.yield_SNII.yield,
+ "feedback-tables", (void **)&feedback_props->yield_SNII.yield,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.SNII_n_z * stars->feedback.SNII_n_mass *
- stars->feedback.SNII_n_elements * sizeof(double)) != 0) {
+ feedback_props->SNII_n_z * feedback_props->SNII_n_mass *
+ feedback_props->SNII_n_elements * sizeof(double)) != 0) {
error("Failed to allocate SNII yield array");
}
/* Allocate array to store SNII yield table resampled by IMF mass bins */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_SNII.yield_IMF_resampled,
+ (void **)&feedback_props->yield_SNII.yield_IMF_resampled,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.SNII_n_z *
- stars->feedback.n_imf_mass_bins *
+ feedback_props->SNII_n_z *
+ feedback_props->n_imf_mass_bins *
chemistry_element_count * sizeof(double)) != 0) {
error("Failed to allocate SNII IMF resampled array");
}
/* Allocate array to store SNII ejecta tables */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_SNII.ejecta,
+ (void **)&feedback_props->yield_SNII.ejecta,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.SNII_n_z * stars->feedback.SNII_n_mass *
+ feedback_props->SNII_n_z * feedback_props->SNII_n_mass *
sizeof(double)) != 0) {
error("Failed to allocate SNII ejecta array");
}
/* Allocate array to store SNII ejecta table resampled by IMF mass bins */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_SNII.ejecta_IMF_resampled,
+ (void **)&feedback_props->yield_SNII.ejecta_IMF_resampled,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.SNII_n_z *
- stars->feedback.n_imf_mass_bins * sizeof(double)) !=
+ feedback_props->SNII_n_z *
+ feedback_props->n_imf_mass_bins * sizeof(double)) !=
0) {
error("Failed to allocate SNII ejecta IMF resampled array");
}
/* Allocate array to store table of total metals released by SNII */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_SNII.total_metals,
+ (void **)&feedback_props->yield_SNII.total_metals,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.SNII_n_z * stars->feedback.SNII_n_mass *
+ feedback_props->SNII_n_z * feedback_props->SNII_n_mass *
sizeof(double)) != 0) {
error("Failed to allocate SNII total metals array");
}
@@ -531,69 +503,69 @@ inline static void allocate_yield_tables(struct stars_props *restrict stars) {
* IMF mass bins */
if (swift_memalign(
"feedback-tables",
- (void **)&stars->feedback.yield_SNII.total_metals_IMF_resampled,
+ (void **)&feedback_props->yield_SNII.total_metals_IMF_resampled,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.SNII_n_z * stars->feedback.n_imf_mass_bins *
+ feedback_props->SNII_n_z * feedback_props->n_imf_mass_bins *
sizeof(double)) != 0) {
error("Failed to allocate SNII total metals IMF resampled array");
}
/* Allocate array for lifetimes mass bins */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.lifetimes.mass,
+ (void **)&feedback_props->lifetimes.mass,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.lifetimes.n_mass * sizeof(double)) != 0) {
+ feedback_props->lifetimes.n_mass * sizeof(double)) != 0) {
error("Failed to allocate lifetime mass array");
}
/* Allocate array for lifetimes metallicity bins */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.lifetimes.metallicity,
+ (void **)&feedback_props->lifetimes.metallicity,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.lifetimes.n_z * sizeof(double)) != 0) {
+ feedback_props->lifetimes.n_z * sizeof(double)) != 0) {
error("Failed to allocate lifetime metallicity array");
}
/* Allocate lifetimes array */
- stars->feedback.lifetimes.dyingtime =
- (double **)malloc(stars->feedback.lifetimes.n_z * sizeof(double *));
- for (int i = 0; i < stars->feedback.lifetimes.n_z; i++) {
- stars->feedback.lifetimes.dyingtime[i] =
- (double *)malloc(stars->feedback.lifetimes.n_mass * sizeof(double));
+ feedback_props->lifetimes.dyingtime =
+ (double **)malloc(feedback_props->lifetimes.n_z * sizeof(double *));
+ for (int i = 0; i < feedback_props->lifetimes.n_z; i++) {
+ feedback_props->lifetimes.dyingtime[i] =
+ (double *)malloc(feedback_props->lifetimes.n_mass * sizeof(double));
}
/* Allocate SNII factor array */
- if (swift_memalign("feedback-tables", (void **)&stars->feedback.typeII_factor,
+ if (swift_memalign("feedback-tables", (void **)&feedback_props->typeII_factor,
SWIFT_STRUCT_ALIGNMENT,
chemistry_element_count * sizeof(double)) != 0) {
error("Failed to allocate SNII factor array");
}
/* Allocate arrays to store names of elements tracked for SNIa, SNII, AGB */
- stars->feedback.SNIa_element_names =
- (char **)malloc(stars->feedback.SNIa_n_elements * sizeof(char *));
- for (int i = 0; i < stars->feedback.SNIa_n_elements; i++) {
- stars->feedback.SNIa_element_names[i] =
- (char *)malloc(stars->feedback.element_name_length * sizeof(char));
+ feedback_props->SNIa_element_names =
+ (char **)malloc(feedback_props->SNIa_n_elements * sizeof(char *));
+ for (int i = 0; i < feedback_props->SNIa_n_elements; i++) {
+ feedback_props->SNIa_element_names[i] =
+ (char *)malloc(feedback_props->element_name_length * sizeof(char));
}
- stars->feedback.SNII_element_names =
- (char **)malloc(stars->feedback.SNII_n_elements * sizeof(char *));
- for (int i = 0; i < stars->feedback.SNII_n_elements; i++) {
- stars->feedback.SNII_element_names[i] =
- (char *)malloc(stars->feedback.element_name_length * sizeof(char));
+ feedback_props->SNII_element_names =
+ (char **)malloc(feedback_props->SNII_n_elements * sizeof(char *));
+ for (int i = 0; i < feedback_props->SNII_n_elements; i++) {
+ feedback_props->SNII_element_names[i] =
+ (char *)malloc(feedback_props->element_name_length * sizeof(char));
}
- stars->feedback.AGB_element_names =
- (char **)malloc(stars->feedback.AGB_n_elements * sizeof(char *));
- for (int i = 0; i < stars->feedback.AGB_n_elements; i++) {
- stars->feedback.AGB_element_names[i] =
- (char *)malloc(stars->feedback.element_name_length * sizeof(char));
+ feedback_props->AGB_element_names =
+ (char **)malloc(feedback_props->AGB_n_elements * sizeof(char *));
+ for (int i = 0; i < feedback_props->AGB_n_elements; i++) {
+ feedback_props->AGB_element_names[i] =
+ (char *)malloc(feedback_props->element_name_length * sizeof(char));
}
/* Allocate array of IMF mass bins */
if (swift_memalign("feedback-tables",
- (void **)&stars->feedback.yield_mass_bins,
+ (void **)&feedback_props->yield_mass_bins,
SWIFT_STRUCT_ALIGNMENT,
- stars->feedback.n_imf_mass_bins * sizeof(double)) != 0) {
+ feedback_props->n_imf_mass_bins * sizeof(double)) != 0) {
error("Failed to allocate imf mass bins array");
}
}
@@ -603,41 +575,41 @@ inline static void allocate_yield_tables(struct stars_props *restrict stars) {
*
* @param stars the #stars_props data structure
*/
-inline static void compute_yields(struct stars_props *restrict stars) {
+inline static void compute_yields(struct feedback_props *feedback_props) {
int flat_index_3d, flat_index_2d;
/* convert SNII tables to log10 */
- for (int i = 0; i < stars->feedback.SNII_n_mass; i++) {
- stars->feedback.yield_SNII.mass[i] =
- log10(stars->feedback.yield_SNII.mass[i]);
- }
- for (int i = 0; i < stars->feedback.SNII_n_z; i++) {
- if (stars->feedback.yield_SNII.metallicity[i] > 0) {
- stars->feedback.yield_SNII.metallicity[i] =
- log10(stars->feedback.yield_SNII.metallicity[i]);
+ for (int i = 0; i < feedback_props->SNII_n_mass; i++) {
+ feedback_props->yield_SNII.mass[i] =
+ log10(feedback_props->yield_SNII.mass[i]);
+ }
+ for (int i = 0; i < feedback_props->SNII_n_z; i++) {
+ if (feedback_props->yield_SNII.metallicity[i] > 0) {
+ feedback_props->yield_SNII.metallicity[i] =
+ log10(feedback_props->yield_SNII.metallicity[i]);
} else {
- stars->feedback.yield_SNII.metallicity[i] = log10_min_metallicity;
+ feedback_props->yield_SNII.metallicity[i] = log10_min_metallicity;
}
}
/* convert AGB tables to log10 */
- for (int i = 0; i < stars->feedback.AGB_n_mass; i++) {
- stars->feedback.yield_AGB.mass[i] =
- log10(stars->feedback.yield_AGB.mass[i]);
- }
- for (int i = 0; i < stars->feedback.AGB_n_z; i++) {
- if (stars->feedback.yield_AGB.metallicity[i] > 0) {
- stars->feedback.yield_AGB.metallicity[i] =
- log10(stars->feedback.yield_AGB.metallicity[i]);
+ for (int i = 0; i < feedback_props->AGB_n_mass; i++) {
+ feedback_props->yield_AGB.mass[i] =
+ log10(feedback_props->yield_AGB.mass[i]);
+ }
+ for (int i = 0; i < feedback_props->AGB_n_z; i++) {
+ if (feedback_props->yield_AGB.metallicity[i] > 0) {
+ feedback_props->yield_AGB.metallicity[i] =
+ log10(feedback_props->yield_AGB.metallicity[i]);
} else {
- stars->feedback.yield_AGB.metallicity[i] = log10_min_metallicity;
+ feedback_props->yield_AGB.metallicity[i] = log10_min_metallicity;
}
}
/* Declare temporary tables to accumulate yields */
- double SNII_yield[stars->feedback.SNII_n_mass];
- double AGB_yield[stars->feedback.AGB_n_mass];
+ double SNII_yield[feedback_props->SNII_n_mass];
+ double AGB_yield[feedback_props->AGB_n_mass];
float result;
/* Resample yields for each element tracked in EAGLE */
@@ -646,101 +618,101 @@ inline static void compute_yields(struct stars_props *restrict stars) {
elem < chemistry_element_count; elem++) {
/* SNIa */
element_index = get_element_index(chemistry_get_element_name(elem),
- stars->feedback.SNIa_element_names,
- stars->feedback.SNIa_n_elements);
- stars->feedback.yield_SNIa_IMF_resampled[elem] =
- stars->feedback.yields_SNIa[element_index];
+ feedback_props->SNIa_element_names,
+ feedback_props->SNIa_n_elements);
+ feedback_props->yield_SNIa_IMF_resampled[elem] =
+ feedback_props->yields_SNIa[element_index];
/* SNII */
element_index = get_element_index(chemistry_get_element_name(elem),
- stars->feedback.SNII_element_names,
- stars->feedback.SNII_n_elements);
- for (int i = 0; i < stars->feedback.SNII_n_z; i++) {
- for (int j = 0; j < stars->feedback.SNII_n_mass; j++) {
+ feedback_props->SNII_element_names,
+ feedback_props->SNII_n_elements);
+ for (int i = 0; i < feedback_props->SNII_n_z; i++) {
+ for (int j = 0; j < feedback_props->SNII_n_mass; j++) {
flat_index_3d = feedback_row_major_index_3d(
- i, element_index, j, stars->feedback.SNII_n_z,
- stars->feedback.SNII_n_elements, stars->feedback.SNII_n_mass);
- SNII_yield[j] = stars->feedback.yield_SNII.yield[flat_index_3d] *
- exp(M_LN10 * (-stars->feedback.yield_SNII.mass[j]));
+ i, element_index, j, feedback_props->SNII_n_z,
+ feedback_props->SNII_n_elements, feedback_props->SNII_n_mass);
+ SNII_yield[j] = feedback_props->yield_SNII.yield[flat_index_3d] *
+ exp(M_LN10 * (-feedback_props->yield_SNII.mass[j]));
}
- for (int k = 0; k < stars->feedback.n_imf_mass_bins; k++) {
- if (stars->feedback.yield_mass_bins[k] <
- stars->feedback.yield_SNII.mass[0])
+ for (int k = 0; k < feedback_props->n_imf_mass_bins; k++) {
+ if (feedback_props->yield_mass_bins[k] <
+ feedback_props->yield_SNII.mass[0])
result = SNII_yield[0];
- else if (stars->feedback.yield_mass_bins[k] >
- stars->feedback.yield_SNII
- .mass[stars->feedback.SNII_n_mass - 1])
- result = SNII_yield[stars->feedback.SNII_n_mass - 1];
+ else if (feedback_props->yield_mass_bins[k] >
+ feedback_props->yield_SNII
+ .mass[feedback_props->SNII_n_mass - 1])
+ result = SNII_yield[feedback_props->SNII_n_mass - 1];
else {
result = interpolate_1D_non_uniform(
- stars->feedback.yield_SNII.mass, SNII_yield,
- stars->feedback.SNII_n_mass, stars->feedback.yield_mass_bins[k]);
+ feedback_props->yield_SNII.mass, SNII_yield,
+ feedback_props->SNII_n_mass, feedback_props->yield_mass_bins[k]);
}
flat_index_3d = feedback_row_major_index_3d(
- i, elem, k, stars->feedback.SNII_n_z, chemistry_element_count,
- stars->feedback.n_imf_mass_bins);
- stars->feedback.yield_SNII.yield_IMF_resampled[flat_index_3d] =
- exp(M_LN10 * stars->feedback.yield_mass_bins[k]) * result;
+ i, elem, k, feedback_props->SNII_n_z, chemistry_element_count,
+ feedback_props->n_imf_mass_bins);
+ feedback_props->yield_SNII.yield_IMF_resampled[flat_index_3d] =
+ exp(M_LN10 * feedback_props->yield_mass_bins[k]) * result;
}
}
- for (int i = 0; i < stars->feedback.SNII_n_z; i++) {
- for (int k = 0; k < stars->feedback.n_imf_mass_bins; k++) {
+ for (int i = 0; i < feedback_props->SNII_n_z; i++) {
+ for (int k = 0; k < feedback_props->n_imf_mass_bins; k++) {
flat_index_2d = feedback_row_major_index_2d(
- i, k, stars->feedback.SNII_n_z, stars->feedback.n_imf_mass_bins);
+ i, k, feedback_props->SNII_n_z, feedback_props->n_imf_mass_bins);
flat_index_3d = feedback_row_major_index_3d(
- i, elem, k, stars->feedback.SNII_n_z, chemistry_element_count,
- stars->feedback.n_imf_mass_bins);
+ i, elem, k, feedback_props->SNII_n_z, chemistry_element_count,
+ feedback_props->n_imf_mass_bins);
if (strcmp(chemistry_get_element_name(elem), "Hydrogen") != 0 ||
strcmp(chemistry_get_element_name(elem), "Helium") != 0) {
- stars->feedback.yield_SNII
+ feedback_props->yield_SNII
.total_metals_IMF_resampled[flat_index_2d] +=
- (stars->feedback.typeII_factor[elem] - 1) *
- stars->feedback.yield_SNII.yield_IMF_resampled[flat_index_3d];
+ (feedback_props->typeII_factor[elem] - 1) *
+ feedback_props->yield_SNII.yield_IMF_resampled[flat_index_3d];
}
- stars->feedback.yield_SNII.yield_IMF_resampled[flat_index_3d] *=
- stars->feedback.typeII_factor[elem];
+ feedback_props->yield_SNII.yield_IMF_resampled[flat_index_3d] *=
+ feedback_props->typeII_factor[elem];
}
}
/* AGB */
element_index = get_element_index(chemistry_get_element_name(elem),
- stars->feedback.AGB_element_names,
- stars->feedback.AGB_n_elements);
+ feedback_props->AGB_element_names,
+ feedback_props->AGB_n_elements);
if (element_index < 0) {
error("element not tracked for AGB");
} else {
- for (int i = 0; i < stars->feedback.AGB_n_z; i++) {
- for (int j = 0; j < stars->feedback.AGB_n_mass; j++) {
+ for (int i = 0; i < feedback_props->AGB_n_z; i++) {
+ for (int j = 0; j < feedback_props->AGB_n_mass; j++) {
flat_index_3d = feedback_row_major_index_3d(
- i, element_index, j, stars->feedback.AGB_n_z,
- stars->feedback.AGB_n_elements, stars->feedback.AGB_n_mass);
- AGB_yield[j] = stars->feedback.yield_AGB.yield[flat_index_3d] *
- exp(M_LN10 * (-stars->feedback.yield_AGB.mass[j]));
+ i, element_index, j, feedback_props->AGB_n_z,
+ feedback_props->AGB_n_elements, feedback_props->AGB_n_mass);
+ AGB_yield[j] = feedback_props->yield_AGB.yield[flat_index_3d] *
+ exp(M_LN10 * (-feedback_props->yield_AGB.mass[j]));
}
- for (int j = 0; j < stars->feedback.n_imf_mass_bins; j++) {
- if (stars->feedback.yield_mass_bins[j] <
- stars->feedback.yield_AGB.mass[0])
+ for (int j = 0; j < feedback_props->n_imf_mass_bins; j++) {
+ if (feedback_props->yield_mass_bins[j] <
+ feedback_props->yield_AGB.mass[0])
result = AGB_yield[0];
- else if (stars->feedback.yield_mass_bins[j] >
- stars->feedback.yield_AGB
- .mass[stars->feedback.AGB_n_mass - 1])
- result = AGB_yield[stars->feedback.AGB_n_mass - 1];
+ else if (feedback_props->yield_mass_bins[j] >
+ feedback_props->yield_AGB
+ .mass[feedback_props->AGB_n_mass - 1])
+ result = AGB_yield[feedback_props->AGB_n_mass - 1];
else
result = interpolate_1D_non_uniform(
- stars->feedback.yield_AGB.mass, AGB_yield,
- stars->feedback.AGB_n_mass, stars->feedback.yield_mass_bins[j]);
+ feedback_props->yield_AGB.mass, AGB_yield,
+ feedback_props->AGB_n_mass, feedback_props->yield_mass_bins[j]);
flat_index_3d = feedback_row_major_index_3d(
- i, elem, j, stars->feedback.AGB_n_z, chemistry_element_count,
- stars->feedback.n_imf_mass_bins);
- stars->feedback.yield_AGB.yield_IMF_resampled[flat_index_3d] =
- exp(M_LN10 * stars->feedback.yield_mass_bins[j]) * result;
+ i, elem, j, feedback_props->AGB_n_z, chemistry_element_count,
+ feedback_props->n_imf_mass_bins);
+ feedback_props->yield_AGB.yield_IMF_resampled[flat_index_3d] =
+ exp(M_LN10 * feedback_props->yield_mass_bins[j]) * result;
}
}
}
@@ -752,123 +724,123 @@ inline static void compute_yields(struct stars_props *restrict stars) {
*
* @param stars the #stars_props data structure
*/
-inline static void compute_ejecta(struct stars_props *restrict stars) {
+inline static void compute_ejecta(struct feedback_props *feedback_props) {
/* Declare temporary tables to accumulate yields */
- double SNII_ejecta[stars->feedback.SNII_n_mass];
- double AGB_ejecta[stars->feedback.AGB_n_mass];
+ double SNII_ejecta[feedback_props->SNII_n_mass];
+ double AGB_ejecta[feedback_props->AGB_n_mass];
float result;
int flat_index;
/* Resample SNII ejecta */
- for (int i = 0; i < stars->feedback.SNII_n_z; i++) {
- for (int k = 0; k < stars->feedback.SNII_n_mass; k++) {
- flat_index = feedback_row_major_index_2d(i, k, stars->feedback.SNII_n_z,
- stars->feedback.SNII_n_mass);
- SNII_ejecta[k] = stars->feedback.yield_SNII.ejecta[flat_index] *
- exp(M_LN10 * (-stars->feedback.yield_SNII.mass[k]));
+ for (int i = 0; i < feedback_props->SNII_n_z; i++) {
+ for (int k = 0; k < feedback_props->SNII_n_mass; k++) {
+ flat_index = feedback_row_major_index_2d(i, k, feedback_props->SNII_n_z,
+ feedback_props->SNII_n_mass);
+ SNII_ejecta[k] = feedback_props->yield_SNII.ejecta[flat_index] *
+ exp(M_LN10 * (-feedback_props->yield_SNII.mass[k]));
}
- for (int k = 0; k < stars->feedback.n_imf_mass_bins; k++) {
- if (stars->feedback.yield_mass_bins[k] <
- stars->feedback.yield_SNII.mass[0])
+ for (int k = 0; k < feedback_props->n_imf_mass_bins; k++) {
+ if (feedback_props->yield_mass_bins[k] <
+ feedback_props->yield_SNII.mass[0])
result = SNII_ejecta[0];
- else if (stars->feedback.yield_mass_bins[k] >
- stars->feedback.yield_SNII.mass[stars->feedback.SNII_n_mass - 1])
- result = SNII_ejecta[stars->feedback.SNII_n_mass - 1];
+ else if (feedback_props->yield_mass_bins[k] >
+ feedback_props->yield_SNII.mass[feedback_props->SNII_n_mass - 1])
+ result = SNII_ejecta[feedback_props->SNII_n_mass - 1];
else
result = interpolate_1D_non_uniform(
- stars->feedback.yield_SNII.mass, SNII_ejecta,
- stars->feedback.SNII_n_mass, stars->feedback.yield_mass_bins[k]);
+ feedback_props->yield_SNII.mass, SNII_ejecta,
+ feedback_props->SNII_n_mass, feedback_props->yield_mass_bins[k]);
- flat_index = feedback_row_major_index_2d(i, k, stars->feedback.SNII_n_z,
- stars->feedback.n_imf_mass_bins);
- stars->feedback.yield_SNII.ejecta_IMF_resampled[flat_index] =
- exp(M_LN10 * stars->feedback.yield_mass_bins[k]) * result;
+ flat_index = feedback_row_major_index_2d(i, k, feedback_props->SNII_n_z,
+ feedback_props->n_imf_mass_bins);
+ feedback_props->yield_SNII.ejecta_IMF_resampled[flat_index] =
+ exp(M_LN10 * feedback_props->yield_mass_bins[k]) * result;
}
}
/* resample SNII total metals released */
- for (int i = 0; i < stars->feedback.SNII_n_z; i++) {
- for (int k = 0; k < stars->feedback.SNII_n_mass; k++) {
- flat_index = feedback_row_major_index_2d(i, k, stars->feedback.SNII_n_z,
- stars->feedback.SNII_n_mass);
- SNII_ejecta[k] = stars->feedback.yield_SNII.total_metals[flat_index] *
- exp(M_LN10 * (-stars->feedback.yield_SNII.mass[k]));
+ for (int i = 0; i < feedback_props->SNII_n_z; i++) {
+ for (int k = 0; k < feedback_props->SNII_n_mass; k++) {
+ flat_index = feedback_row_major_index_2d(i, k, feedback_props->SNII_n_z,
+ feedback_props->SNII_n_mass);
+ SNII_ejecta[k] = feedback_props->yield_SNII.total_metals[flat_index] *
+ exp(M_LN10 * (-feedback_props->yield_SNII.mass[k]));
}
- for (int k = 0; k < stars->feedback.n_imf_mass_bins; k++) {
- if (stars->feedback.yield_mass_bins[k] <
- stars->feedback.yield_SNII.mass[0])
+ for (int k = 0; k < feedback_props->n_imf_mass_bins; k++) {
+ if (feedback_props->yield_mass_bins[k] <
+ feedback_props->yield_SNII.mass[0])
result = SNII_ejecta[0];
- else if (stars->feedback.yield_mass_bins[k] >
- stars->feedback.yield_SNII.mass[stars->feedback.SNII_n_mass - 1])
- result = SNII_ejecta[stars->feedback.SNII_n_mass - 1];
+ else if (feedback_props->yield_mass_bins[k] >
+ feedback_props->yield_SNII.mass[feedback_props->SNII_n_mass - 1])
+ result = SNII_ejecta[feedback_props->SNII_n_mass - 1];
else
result = interpolate_1D_non_uniform(
- stars->feedback.yield_SNII.mass, SNII_ejecta,
- stars->feedback.SNII_n_mass, stars->feedback.yield_mass_bins[k]);
+ feedback_props->yield_SNII.mass, SNII_ejecta,
+ feedback_props->SNII_n_mass, feedback_props->yield_mass_bins[k]);
- flat_index = feedback_row_major_index_2d(i, k, stars->feedback.SNII_n_z,
- stars->feedback.n_imf_mass_bins);
- stars->feedback.yield_SNII.total_metals_IMF_resampled[flat_index] =
- exp(M_LN10 * stars->feedback.yield_mass_bins[k]) * result;
+ flat_index = feedback_row_major_index_2d(i, k, feedback_props->SNII_n_z,
+ feedback_props->n_imf_mass_bins);
+ feedback_props->yield_SNII.total_metals_IMF_resampled[flat_index] =
+ exp(M_LN10 * feedback_props->yield_mass_bins[k]) * result;
}
}
/* AGB yields */
- for (int i = 0; i < stars->feedback.AGB_n_z; i++) {
- for (int k = 0; k < stars->feedback.AGB_n_mass; k++) {
- flat_index = feedback_row_major_index_2d(i, k, stars->feedback.AGB_n_z,
- stars->feedback.AGB_n_mass);
- AGB_ejecta[k] = stars->feedback.yield_AGB.ejecta[flat_index] /
- exp(M_LN10 * stars->feedback.yield_AGB.mass[k]);
+ for (int i = 0; i < feedback_props->AGB_n_z; i++) {
+ for (int k = 0; k < feedback_props->AGB_n_mass; k++) {
+ flat_index = feedback_row_major_index_2d(i, k, feedback_props->AGB_n_z,
+ feedback_props->AGB_n_mass);
+ AGB_ejecta[k] = feedback_props->yield_AGB.ejecta[flat_index] /
+ exp(M_LN10 * feedback_props->yield_AGB.mass[k]);
}
- for (int k = 0; k < stars->feedback.n_imf_mass_bins; k++) {
- if (stars->feedback.yield_mass_bins[k] <
- stars->feedback.yield_AGB.mass[0])
+ for (int k = 0; k < feedback_props->n_imf_mass_bins; k++) {
+ if (feedback_props->yield_mass_bins[k] <
+ feedback_props->yield_AGB.mass[0])
result = AGB_ejecta[0];
- else if (stars->feedback.yield_mass_bins[k] >
- stars->feedback.yield_AGB.mass[stars->feedback.AGB_n_mass - 1])
- result = AGB_ejecta[stars->feedback.AGB_n_mass - 1];
+ else if (feedback_props->yield_mass_bins[k] >
+ feedback_props->yield_AGB.mass[feedback_props->AGB_n_mass - 1])
+ result = AGB_ejecta[feedback_props->AGB_n_mass - 1];
else
result = interpolate_1D_non_uniform(
- stars->feedback.yield_AGB.mass, AGB_ejecta,
- stars->feedback.AGB_n_mass, stars->feedback.yield_mass_bins[k]);
+ feedback_props->yield_AGB.mass, AGB_ejecta,
+ feedback_props->AGB_n_mass, feedback_props->yield_mass_bins[k]);
- flat_index = feedback_row_major_index_2d(i, k, stars->feedback.AGB_n_z,
- stars->feedback.n_imf_mass_bins);
- stars->feedback.yield_AGB.ejecta_IMF_resampled[flat_index] =
- exp(M_LN10 * stars->feedback.yield_mass_bins[k]) * result;
+ flat_index = feedback_row_major_index_2d(i, k, feedback_props->AGB_n_z,
+ feedback_props->n_imf_mass_bins);
+ feedback_props->yield_AGB.ejecta_IMF_resampled[flat_index] =
+ exp(M_LN10 * feedback_props->yield_mass_bins[k]) * result;
}
}
- for (int i = 0; i < stars->feedback.AGB_n_z; i++) {
- for (int k = 0; k < stars->feedback.AGB_n_mass; k++) {
- flat_index = feedback_row_major_index_2d(i, k, stars->feedback.AGB_n_z,
- stars->feedback.AGB_n_mass);
- AGB_ejecta[k] = stars->feedback.yield_AGB.total_metals[flat_index] *
- exp(M_LN10 * (-stars->feedback.yield_AGB.mass[k]));
+ for (int i = 0; i < feedback_props->AGB_n_z; i++) {
+ for (int k = 0; k < feedback_props->AGB_n_mass; k++) {
+ flat_index = feedback_row_major_index_2d(i, k, feedback_props->AGB_n_z,
+ feedback_props->AGB_n_mass);
+ AGB_ejecta[k] = feedback_props->yield_AGB.total_metals[flat_index] *
+ exp(M_LN10 * (-feedback_props->yield_AGB.mass[k]));
}
- for (int k = 0; k < stars->feedback.n_imf_mass_bins; k++) {
- if (stars->feedback.yield_mass_bins[k] <
- stars->feedback.yield_AGB.mass[0])
+ for (int k = 0; k < feedback_props->n_imf_mass_bins; k++) {
+ if (feedback_props->yield_mass_bins[k] <
+ feedback_props->yield_AGB.mass[0])
result = AGB_ejecta[0];
- else if (stars->feedback.yield_mass_bins[k] >
- stars->feedback.yield_AGB.mass[stars->feedback.AGB_n_mass - 1])
- result = AGB_ejecta[stars->feedback.AGB_n_mass - 1];
+ else if (feedback_props->yield_mass_bins[k] >
+ feedback_props->yield_AGB.mass[feedback_props->AGB_n_mass - 1])
+ result = AGB_ejecta[feedback_props->AGB_n_mass - 1];
else
result = interpolate_1D_non_uniform(
- stars->feedback.yield_AGB.mass, AGB_ejecta,
- stars->feedback.AGB_n_mass, stars->feedback.yield_mass_bins[k]);
+ feedback_props->yield_AGB.mass, AGB_ejecta,
+ feedback_props->AGB_n_mass, feedback_props->yield_mass_bins[k]);
- flat_index = feedback_row_major_index_2d(i, k, stars->feedback.AGB_n_z,
- stars->feedback.n_imf_mass_bins);
- stars->feedback.yield_AGB.total_metals_IMF_resampled[flat_index] =
- exp(M_LN10 * stars->feedback.yield_mass_bins[k]) * result;
+ flat_index = feedback_row_major_index_2d(i, k, feedback_props->AGB_n_z,
+ feedback_props->n_imf_mass_bins);
+ feedback_props->yield_AGB.total_metals_IMF_resampled[flat_index] =
+ exp(M_LN10 * feedback_props->yield_mass_bins[k]) * result;
}
}
}
diff --git a/src/feedback_struct.h b/src/feedback_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..be6d480e7a7a6165fd6f73e4ca0ed15078f06488
--- /dev/null
+++ b/src/feedback_struct.h
@@ -0,0 +1,39 @@
+/*******************************************************************************
+ * 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 <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_FEEDBACK_STRUCT_H
+#define SWIFT_FEEDBACK_STRUCT_H
+
+/**
+ * @file src/feedback_struct.h
+ * @brief Branches between the different feedback functions.
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Import the right feedback definition */
+#if defined(FEEDBACK_NONE)
+#include "./feedback/none/feedback_struct.h"
+#elif defined(FEEDBACK_EAGLE)
+#include "./feedback/EAGLE/feedback_struct.h"
+#else
+#error "Invalid choice of feedback function."
+#endif
+
+#endif /* SWIFT_FEEDBACK_STRUCT_H */
diff --git a/src/runner.c b/src/runner.c
index 30ae3aa64cbffac6c16a1c0a75eca8f50a79d913..7959a6fe736679bed2b4c5da685d0fbe07dba152 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -50,6 +50,7 @@
#include "engine.h"
#include "entropy_floor.h"
#include "error.h"
+#include "feedback.h"
#include "gravity.h"
#include "hydro.h"
#include "hydro_properties.h"
diff --git a/src/stars.h b/src/stars.h
index bef3cecbd11f26dcd6260aa1f75f3fd37013e96d..dd8390e0206580fc2a07a08e51bb69c6ee5ab5ed 100644
--- a/src/stars.h
+++ b/src/stars.h
@@ -23,20 +23,14 @@
#include "../config.h"
/* Select the correct star model */
-#if defined(FEEDBACK_CONST)
-#include "./stars/const/stars.h"
-#include "./stars/const/stars_iact.h"
-#elif defined(STARS_NONE)
+#if defined(STARS_NONE)
#include "./stars/Default/stars.h"
#include "./stars/Default/stars_iact.h"
#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
-#endif
+#endif /* SWIFT_STARS_H */
diff --git a/src/stars/EAGLE/stars_part.h b/src/stars/EAGLE/stars_part.h
index c9dafee7e1e897f4b83d34c8326dcd1dcbd14a23..7bf82d03ecd056e08a339fcca59f891ca74b7f80 100644
--- a/src/stars/EAGLE/stars_part.h
+++ b/src/stars/EAGLE/stars_part.h
@@ -25,6 +25,7 @@
/* Read chemistry */
#include "chemistry_struct.h"
+#include "feedback_struct.h"
#include "tracers_struct.h"
/**
@@ -97,6 +98,9 @@ struct spart {
/*! Birth density */
float birth_density;
+ /*! Feedback structure */
+ struct feedback_part_data feedback_data;
+
/*! Tracer structure */
struct tracers_xpart_data tracers_data;