diff --git a/configure.ac b/configure.ac
index a43dc8bd915770fa0d46422e755bae43fa0f6f07..2a15c2190d68dfda1d517a98d77057d0741a669f 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1760,9 +1760,6 @@ case "$with_feedback" in
thermal)
AC_DEFINE([FEEDBACK_THERMAL], [1], [Thermal Blastwave])
;;
- const)
- AC_DEFINE([FEEDBACK_CONST], [1], [Constant feedback model])
- ;;
none)
;;
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index ada36234d6a9691a8d78a33a047088920e2c6a11..3e100435df8e75f08b2baa46f2736b015b619147 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -105,3 +105,10 @@ EAGLEEntropyFloor:
Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
+
+EagleStellarEvolution:
+ filename: /cosma5/data/Eagle/BG_Tables/YieldTables/
+ imf_model: Chabrier
+
+EAGLEFeedback:
+ lifetime_flag: 2
diff --git a/src/Makefile.am b/src/Makefile.am
index a7dfecde337b2e87adb52aa70ef5d339e1187cf1..d7e4249a7ff67132505e3a7df8a134d4cd8b266c 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -138,8 +138,6 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
stars.h stars_io.h \
stars/Default/stars.h stars/Default/stars_iact.h stars/Default/stars_io.h \
stars/Default/stars_debug.h stars/Default/stars_part.h \
- stars/const/stars.h stars/const/stars_iact.h stars/const/stars_io.h \
- stars/const/stars_debug.h stars/const/stars_part.h \
stars/EAGLE/stars.h stars/EAGLE/stars_iact.h stars/EAGLE/stars_io.h \
stars/EAGLE/stars_debug.h stars/EAGLE/stars_part.h \
potential/none/potential.h potential/point_mass/potential.h \
diff --git a/src/stars/const/stars.h b/src/stars/const/stars.h
deleted file mode 100644
index 28debcdaf4852b76c96eaad0c08f8a09cfbcc27f..0000000000000000000000000000000000000000
--- a/src/stars/const/stars.h
+++ /dev/null
@@ -1,256 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Coypright (c) 2016 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_CONST_STARS_H
-#define SWIFT_CONST_STARS_H
-
-#include <float.h>
-#include "minmax.h"
-
-/**
- * @brief Computes the gravity time-step of a given star particle.
- *
- * @param sp Pointer to the s-particle data.
- */
-__attribute__((always_inline)) INLINE static float stars_compute_timestep(
- const struct spart* const sp) {
-
- return FLT_MAX;
-}
-
-/**
- * @brief Initialises the s-particles for the first time
- *
- * This function is called only once just after the ICs have been
- * read in to do some conversions.
- *
- * @param sp The particle to act upon
- */
-__attribute__((always_inline)) INLINE static void stars_first_init_spart(
- struct spart* sp) {
-
- sp->time_bin = 0;
-}
-
-/**
- * @brief Prepares a s-particle for its interactions
- *
- * @param sp The particle to act upon
- */
-__attribute__((always_inline)) INLINE static void stars_init_spart(
- struct spart* sp) {
-
-#ifdef DEBUG_INTERACTIONS_STARS
- for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
- sp->ids_ngbs_density[i] = -1;
- sp->num_ngb_density = 0;
-#endif
-
- sp->density.wcount = 0.f;
- sp->density.wcount_dh = 0.f;
-
- sp->omega_normalisation_inv = 0.f;
- sp->ngb_mass = 0.f;
-}
-
-/**
- * @brief Predict additional particle fields forward in time when drifting
- *
- * @param p The particle
- * @param dt_drift The drift time-step for positions.
- */
-__attribute__((always_inline)) INLINE static void stars_predict_extra(
- struct spart* restrict sp, float dt_drift) {
-
- const float h_inv = 1.f / sp->h;
-
- /* Predict smoothing length */
- const float w1 = sp->feedback.h_dt * h_inv * dt_drift;
- if (fabsf(w1) < 0.2f)
- sp->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
- else
- sp->h *= expf(w1);
-}
-
-/**
- * @brief Sets the values to be predicted in the drifts to their values at a
- * kick time
- *
- * @param sp The particle.
- */
-__attribute__((always_inline)) INLINE static void stars_reset_predicted_values(
- struct spart* restrict sp) {}
-
-/**
- * @brief Finishes the calculation of feedback acting on stars
- *
- * Multiplies the forces and accelerations by the appropiate constants
- *
- * @param sp The particle to act upon
- */
-__attribute__((always_inline)) INLINE static void stars_end_feedback(
- struct spart* sp) {}
-
-/**
- * @brief Kick the additional variables
- *
- * @param sp The particle to act upon
- * @param dt The time-step for this kick
- */
-__attribute__((always_inline)) INLINE static void stars_kick_extra(
- struct spart* sp, float dt) {}
-
-/**
- * @brief Finishes the calculation of density on stars
- *
- * @param sp The particle to act upon
- * @param cosmo The current cosmological model.
- */
-__attribute__((always_inline)) INLINE static void stars_end_density(
- struct spart* sp, const struct cosmology* cosmo) {
-
- /* Some smoothing length multiples. */
- const float h = sp->h;
- const float h_inv = 1.0f / h; /* 1/h */
- const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
- const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
-
- /* Finish the calculation by inserting the missing h-factors */
- sp->density.wcount *= h_inv_dim;
- sp->density.wcount_dh *= h_inv_dim_plus_one;
-}
-
-/**
- * @brief Sets all particle fields to sensible values when the #spart has 0
- * ngbs.
- *
- * @param sp The particle to act upon
- * @param cosmo The current cosmological model.
- */
-__attribute__((always_inline)) INLINE static void stars_spart_has_no_neighbours(
- struct spart* restrict sp, const struct cosmology* cosmo) {
-
- /* Some smoothing length multiples. */
- const float h = sp->h;
- const float h_inv = 1.0f / h; /* 1/h */
- const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
-
- /* Re-set problematic values */
- sp->density.wcount = kernel_root * h_inv_dim;
- sp->density.wcount_dh = 0.f;
-}
-
-/**
- * @brief Reset acceleration fields of a particle
- *
- * This is the equivalent of hydro_reset_acceleration.
- * We do not compute the acceleration on star, therefore no need to use it.
- *
- * @param p The particle to act upon
- */
-__attribute__((always_inline)) INLINE static void stars_reset_acceleration(
- struct spart* restrict p) {
-#ifdef DEBUG_INTERACTIONS_STARS
- p->num_ngb_force = 0;
-#endif
-}
-
-/**
- * @brief Evolve the stellar properties of a #spart.
- *
- * This function allows for example to compute the SN rate before sending
- * this information to a different MPI rank.
- *
- * @param sp The particle to act upon
- * @param cosmo The current cosmological model.
- * @param stars_properties The #stars_props
- * @param dt Timestep over which the particle evolves.
- */
-__attribute__((always_inline)) INLINE static void stars_evolve_spart(
- struct spart* restrict sp, const struct stars_props* stars_properties,
- const struct cosmology* cosmo, const struct unit_system* us,
- float current_time, double dt) {
-
- /* Proportion of quantities to be released each timestep */
- float feedback_factor = dt / stars_properties->feedback_timescale;
-
- /* Amount of mass to distribute in one step */
- sp->to_distribute.mass = sp->mass_init * feedback_factor;
-
- /* Set all enrichment quantities to constant values */
- for (int i = 0; i < chemistry_element_count; i++)
- sp->to_distribute.chemistry_data.metal_mass_fraction[i] =
- 1.f / chemistry_element_count;
- sp->to_distribute.chemistry_data.metal_mass_fraction_total =
- 1.f - 2.f / chemistry_element_count;
- sp->to_distribute.chemistry_data.mass_from_AGB =
- 1.0e-2 * sp->to_distribute.mass;
- sp->to_distribute.chemistry_data.metal_mass_fraction_from_AGB = 1.0e-2;
- sp->to_distribute.chemistry_data.mass_from_SNII =
- 1.0e-2 * sp->to_distribute.mass;
- sp->to_distribute.chemistry_data.metal_mass_fraction_from_SNII = 1.0e-2;
- sp->to_distribute.chemistry_data.mass_from_SNIa =
- 1.0e-2 * sp->to_distribute.mass;
- sp->to_distribute.chemistry_data.metal_mass_fraction_from_SNIa = 1.0e-2;
- sp->to_distribute.chemistry_data.iron_mass_fraction_from_SNIa = 1.0e-2;
-
- /* Set feedback to constant values */
- const float total_sn = sp->mass_init / stars_properties->const_solar_mass *
- stars_properties->sn_per_msun;
-
- /* Print total_sn and timescale to be read by test script for checking
- * stochastic energy injection */
- if (dt == 0) {
- FILE* feedback_output = fopen("feedback_properties.dat", "w");
- fprintf(feedback_output, "%.5e \n%.5e \n", total_sn,
- stars_properties->feedback_timescale);
- fclose(feedback_output);
- }
-
- /* Set number of supernovae to distribute */
- sp->to_distribute.num_SNIa = total_sn * feedback_factor;
-
- /* Set ejected thermal energy */
- sp->to_distribute.ejecta_specific_thermal_energy = 1.0e-3;
-}
-
-inline static void stars_evolve_init(struct swift_params* params,
- struct stars_props* restrict stars) {}
-
-/**
- * @brief Reset acceleration fields of a particle
- *
- * This is the equivalent of hydro_reset_acceleration.
- * We do not compute the acceleration on star, therefore no need to use it.
- *
- * @param p The particle to act upon
- */
-__attribute__((always_inline)) INLINE static void stars_reset_feedback(
- struct spart* restrict p) {
-
- /* Reset time derivative */
- p->feedback.h_dt = 0.f;
-
-#ifdef DEBUG_INTERACTIONS_STARS
- for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
- p->ids_ngbs_force[i] = -1;
- p->num_ngb_force = 0;
-#endif
-}
-
-#endif /* SWIFT_CONST_STARS_H */
diff --git a/src/stars/const/stars_debug.h b/src/stars/const/stars_debug.h
deleted file mode 100644
index 53a54b028ad92fce26a41ba46569db745a51544d..0000000000000000000000000000000000000000
--- a/src/stars/const/stars_debug.h
+++ /dev/null
@@ -1,31 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Coypright (c) 2016 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_CONST_STARS_DEBUG_H
-#define SWIFT_CONST_STARS_DEBUG_H
-
-__attribute__((always_inline)) INLINE static void stars_debug_particle(
- const struct spart* p) {
- printf(
- "x=[%.3e,%.3e,%.3e], "
- "v_full=[%.3e,%.3e,%.3e] p->mass=%.3e \n t_begin=%d, t_end=%d\n",
- p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2],
- p->mass, p->ti_begin, p->ti_end);
-}
-
-#endif /* SWIFT_CONST_STARS_DEBUG_H */
diff --git a/src/stars/const/stars_iact.h b/src/stars/const/stars_iact.h
deleted file mode 100644
index c0a2fb3f429043fe3134e6c671f5d5594a63744e..0000000000000000000000000000000000000000
--- a/src/stars/const/stars_iact.h
+++ /dev/null
@@ -1,262 +0,0 @@
-#include "random.h"
-
-/**
- * @brief Density interaction between two particles (non-symmetric).
- *
- * @param r2 Comoving square distance between the two particles.
- * @param dx Comoving vector separating both particles (pi - pj).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param si First sparticle.
- * @param pj Second particle (not updated).
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- * @param xp Extra particle data
- * @param ti_current Current integer time value
- */
-__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_stars_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,
- struct xpart *restrict xp, integertime_t ti_current) {
-
- float wi, wi_dx;
-
- /* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
-
- /* Compute the kernel function */
- const float hi_inv = 1.0f / hi;
- const float ui = r * hi_inv;
- kernel_deval(ui, &wi, &wi_dx);
-
- float wj, wj_dx;
- const float hj_inv = 1.0f / hj;
- const float uj = r * hj_inv;
- kernel_deval(uj, &wj, &wj_dx);
-
- /* Compute contribution to the number of neighbours */
- si->density.wcount += wi;
- si->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
-
- /* Add mass of pj to neighbour mass of si */
- si->ngb_mass += hydro_get_mass(pj);
-
- /* Add contribution of pj to normalisation of kernel (TODO: IMPROVE COMMENT?)
- */
- si->omega_normalisation_inv += wj / hydro_get_physical_density(pj, cosmo);
-
-#ifdef DEBUG_INTERACTIONS_STARS
- /* Update ngb counters */
- if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_STARS)
- si->ids_ngbs_density[si->num_ngb_density] = pj->id;
- ++si->num_ngb_density;
-#endif
-}
-
-/**
- * @brief Increases thermal energy of particle due
- * to feedback by specified amount
- *
- * @param du change in internal energy
- * @param p Particle we're acting on
- * @param xp Extra particle data
- * @param cosmo Cosmology struct
- */
-static inline void thermal_feedback(float du, struct part *restrict p,
- struct xpart *restrict xp,
- const struct cosmology *restrict cosmo) {
-
- float u = hydro_get_physical_internal_energy(p, xp, cosmo);
- hydro_set_physical_internal_energy(p, cosmo, u + du);
- // Just setting p->entropy is not enough because xp->entropy_full gets updated
- // with p->entropy_dt
- // TODO: ADD HYDRO FUNCTIONS FOR UPDATING DRIFTED AND NON DRIFTED INTERNAL
- // ENERGY AND GET RID OF THE ENTROPY UPDATE HERE.
- xp->entropy_full = p->entropy;
-}
-
-/**
- * @brief Feedback interaction between two particles (non-symmetric).
- * Used for updating properties of gas particles neighbouring a star particle
- *
- * @param r2 Comoving square distance between the two particles.
- * @param dx Comoving vector separating both particles (si - pj).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param si First (star) particle.
- * @param pj Second (gas) particle.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- * @param xp Extra particle data
- * @param ti_current Current integer time used value for seeding random number
- * generator
- */
-__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_stars_feedback(
- float r2, const float *dx, float hi, float hj, struct spart *restrict si,
- struct part *restrict pj, const struct cosmology *restrict cosmo,
- const struct stars_props *restrict stars_properties,
- struct xpart *restrict xp, integertime_t ti_current) {
- float wj;
-
- /* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
-
- /* Compute the kernel function */
- const float hj_inv = 1.0f / hj;
- const float uj = r * hj_inv;
- kernel_eval(uj, &wj);
-
- /* Compute weighting for distributing feedback quantities */
- const float omega_frac = wj / (hydro_get_physical_density(pj, cosmo) *
- si->omega_normalisation_inv);
-
- /* Update particle mass */
- const float current_mass = hydro_get_mass(pj);
- float new_mass = current_mass + si->to_distribute.mass * omega_frac;
- if (stars_properties->const_feedback_energy_testing) new_mass = current_mass;
- hydro_set_mass(pj, new_mass);
-
- message(
- "particle %llu current mass %.5e new mass %.5e mass to distribute %.5e "
- "omega_frac %.5e norm %.5e ngb_mass %.5e",
- pj->id, current_mass, new_mass, si->to_distribute.mass, omega_frac,
- si->omega_normalisation_inv, si->ngb_mass);
-
- /* Decrease the mass of star particle */
- si->mass -= si->to_distribute.mass * omega_frac;
-
- // ALEXEI: do we want to use the smoothed mass fraction?
- /* Update total metallicity */
- const float current_metal_mass_total =
- pj->chemistry_data.metal_mass_fraction_total * current_mass;
- const float new_metal_mass_total =
- current_metal_mass_total +
- si->to_distribute.chemistry_data.metal_mass_fraction_total *
- si->to_distribute.mass * omega_frac;
- pj->chemistry_data.metal_mass_fraction_total =
- new_metal_mass_total / new_mass;
-
- /* Update mass fraction of each tracked element */
- for (int elem = 0; elem < chemistry_element_count; elem++) {
- const float current_metal_mass =
- pj->chemistry_data.metal_mass_fraction[elem] * current_mass;
- const float new_metal_mass =
- current_metal_mass +
- si->to_distribute.chemistry_data.metal_mass_fraction[elem] *
- si->to_distribute.mass * omega_frac;
- pj->chemistry_data.metal_mass_fraction[elem] = new_metal_mass / new_mass;
- }
-
- /* Update iron mass fraction from SNIa */
- const float current_iron_from_SNIa_mass =
- pj->chemistry_data.iron_mass_fraction_from_SNIa * current_mass;
- const float new_iron_from_SNIa_mass =
- current_iron_from_SNIa_mass +
- si->to_distribute.chemistry_data.iron_mass_fraction_from_SNIa *
- si->to_distribute.mass * omega_frac;
- pj->chemistry_data.iron_mass_fraction_from_SNIa =
- new_iron_from_SNIa_mass / new_mass;
-
- /* Update mass from SNIa */
- pj->chemistry_data.mass_from_SNIa +=
- si->to_distribute.chemistry_data.mass_from_SNIa * omega_frac;
-
- /* Update metal mass fraction from SNIa */
- const float current_metal_mass_fraction_from_SNIa =
- pj->chemistry_data.metal_mass_fraction_from_SNIa * current_mass;
- const float new_metal_mass_fraction_from_SNIa =
- current_metal_mass_fraction_from_SNIa +
- si->to_distribute.chemistry_data.metal_mass_fraction_from_SNIa *
- si->to_distribute.mass * omega_frac;
- pj->chemistry_data.metal_mass_fraction_from_SNIa =
- new_metal_mass_fraction_from_SNIa / new_mass;
-
- /* Update mass from SNII */
- pj->chemistry_data.mass_from_SNII +=
- si->to_distribute.chemistry_data.mass_from_SNII * omega_frac;
-
- /* Update metal mass fraction from SNII */
- const float current_metal_mass_fraction_from_SNII =
- pj->chemistry_data.metal_mass_fraction_from_SNII * current_mass;
- const float new_metal_mass_fraction_from_SNII =
- current_metal_mass_fraction_from_SNII +
- si->to_distribute.chemistry_data.metal_mass_fraction_from_SNII *
- si->to_distribute.mass * omega_frac;
- pj->chemistry_data.metal_mass_fraction_from_SNII =
- new_metal_mass_fraction_from_SNII / new_mass;
-
- /* Update mass from AGB */
- pj->chemistry_data.mass_from_AGB +=
- si->to_distribute.chemistry_data.mass_from_AGB * omega_frac;
-
- /* Update metal mass fraction from AGB */
- const float current_metal_mass_fraction_from_AGB =
- pj->chemistry_data.metal_mass_fraction_from_AGB * current_mass;
- const float new_metal_mass_fraction_from_AGB =
- current_metal_mass_fraction_from_AGB +
- si->to_distribute.chemistry_data.metal_mass_fraction_from_AGB *
- si->to_distribute.mass * omega_frac;
- pj->chemistry_data.metal_mass_fraction_from_AGB =
- new_metal_mass_fraction_from_AGB / new_mass;
-
- /* Update momentum */
- for (int i = 0; i < 3; i++) {
- // Do we need to calculate relative velocities here?
- pj->v[i] += si->to_distribute.mass * omega_frac * si->v[i];
- }
-
- /* Energy feedback */
- float heating_probability = -1.f, du = 0.f, d_energy = 0.f;
- d_energy =
- si->to_distribute.mass *
- (si->to_distribute.ejecta_specific_thermal_energy +
- 0.5 * (si->v[0] * si->v[0] + si->v[1] * si->v[1] + si->v[2] * si->v[2]) *
- cosmo->a2_inv);
-
- // If statement temporary for testing, in practice would always be on.
- if (stars_properties->const_feedback_energy_testing) {
- if (stars_properties->continuous_heating) {
- // We're doing ONLY continuous heating
- d_energy += si->to_distribute.num_SNIa *
- stars_properties->total_energy_SNe * omega_frac;
- du = d_energy / hydro_get_mass(pj);
- if (du > 0)
- message("id %llu du %.5e initial u %.5e", pj->id, du,
- hydro_get_physical_internal_energy(pj, xp, cosmo));
- thermal_feedback(du, pj, xp, cosmo);
- } else {
- // We're doing stochastic heating
- heating_probability =
- stars_properties->SNe_temperature * si->to_distribute.num_SNIa *
- stars_properties->SNIa_energy_fraction /
- (stars_properties->SNe_deltaT_desired * si->ngb_mass);
- du = stars_properties->SNe_deltaT_desired *
- stars_properties->temp_to_u_factor;
- if (heating_probability >= 1) {
- du = stars_properties->total_energy_SNe * si->to_distribute.num_SNIa /
- si->ngb_mass;
- heating_probability = 1;
- }
- }
- }
-
- double random_num =
- random_unit_interval(pj->id, ti_current, random_number_stellar_feedback);
- if (random_num < heating_probability) {
- message(
- "we did some heating! id %llu probability %.5e random_num %.5e du %.5e "
- "initial u %.5e",
- pj->id, heating_probability, random_num, du,
- hydro_get_physical_internal_energy(pj, xp, cosmo));
- thermal_feedback(du, pj, xp, cosmo);
- FILE *feedback_output = fopen("feedback_properties.dat", "a");
- fprintf(feedback_output, "%.5e %.5e\n", heating_probability,
- du * hydro_get_mass(pj));
- fclose(feedback_output);
- }
-}
diff --git a/src/stars/const/stars_io.h b/src/stars/const/stars_io.h
deleted file mode 100644
index d3f184f9aac6592f5b40b4b9fc917020f9eac6e3..0000000000000000000000000000000000000000
--- a/src/stars/const/stars_io.h
+++ /dev/null
@@ -1,253 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Coypright (c) 2016 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_CONST_STARS_IO_H
-#define SWIFT_CONST_STARS_IO_H
-
-#include "io_properties.h"
-#include "stars_part.h"
-
-/**
- * @brief Specifies which s-particle fields to read from a dataset
- *
- * @param sparts The s-particle array.
- * @param list The list of i/o properties to read.
- * @param num_fields The number of i/o fields to read.
- */
-INLINE static void stars_read_particles(struct spart *sparts,
- struct io_props *list,
- int *num_fields) {
-
- /* Say how much we want to read */
- *num_fields = 6;
-
- /* List what we want to read */
- list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
- UNIT_CONV_LENGTH, sparts, x);
- list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
- UNIT_CONV_SPEED, sparts, v);
- list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
- sparts, mass);
- /* Temporary addition, discuss with Folkert what to do about initial mass
- * when sparts are read from ICs. */
- list[3] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
- sparts, mass_init);
- list[4] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY,
- UNIT_CONV_NO_UNITS, sparts, id);
- list[5] = io_make_input_field("SmoothingLength", FLOAT, 1, OPTIONAL,
- UNIT_CONV_LENGTH, sparts, h);
-}
-
-/**
- * @brief Specifies which s-particle fields to write to a dataset
- *
- * @param sparts The s-particle array.
- * @param list The list of i/o properties to write.
- * @param num_fields The number of i/o fields to write.
- */
-INLINE static void stars_write_particles(const struct spart *sparts,
- struct io_props *list,
- int *num_fields) {
-
- /* Say how much we want to read */
- *num_fields = 5;
-
- /* List what we want to read */
- list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
- sparts, x);
- list[1] =
- io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, sparts, v);
- list[2] =
- io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, sparts, mass);
- list[3] = io_make_output_field("ParticleIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS,
- sparts, id);
- list[4] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
- sparts, h);
-}
-
-/**
- * @brief Initialize the global properties of the stars scheme.
- *
- * By default, takes the values provided by the hydro.
- *
- * @param sp The #stars_props.
- * @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.
- * @param cosmo The already read-in cosmology properties.
- */
-INLINE static void stars_props_init(struct stars_props *sp,
- const struct phys_const *phys_const,
- const struct unit_system *us,
- struct swift_params *params,
- const struct hydro_props *p,
- const struct cosmology *cosmo) {
-
- /* Kernel properties */
- sp->eta_neighbours = parser_get_opt_param_float(
- params, "Stars:resolution_eta", p->eta_neighbours);
-
- /* Tolerance for the smoothing length Newton-Raphson scheme */
- sp->h_tolerance =
- parser_get_opt_param_float(params, "Stars:h_tolerance", p->h_tolerance);
-
- /* Get derived properties */
- sp->target_neighbours = pow_dimension(sp->eta_neighbours) * kernel_norm;
- const float delta_eta = sp->eta_neighbours * (1.f + sp->h_tolerance);
- sp->delta_neighbours =
- (pow_dimension(delta_eta) - pow_dimension(sp->eta_neighbours)) *
- kernel_norm;
-
- /* Maximal smoothing length */
- sp->h_max = parser_get_opt_param_float(params, "Stars:h_max", p->h_max);
-
- /* Number of iterations to converge h */
- sp->max_smoothing_iterations = parser_get_opt_param_int(
- params, "Stars:max_ghost_iterations", p->max_smoothing_iterations);
-
- /* Time integration properties */
- const float max_volume_change =
- parser_get_opt_param_float(params, "Stars:max_volume_change", -1);
- if (max_volume_change == -1)
- sp->log_max_h_change = p->log_max_h_change;
- else
- sp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
-
- /* Check if we are heating continuously. Set to 1 if using continuous, 0 for
- * stochastic */
- sp->continuous_heating =
- parser_get_opt_param_int(params, "Stars:continuous_heating", 0);
-
- /* Are we testing the energy injection in the constant feedback model? */
- sp->const_feedback_energy_testing =
- parser_get_opt_param_int(params, "Stars:energy_testing", 0);
-
- /* Set temperature increase due to supernovae */
- sp->SNe_deltaT_desired =
- 3.16228e7 / units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
-
- /* Calculate temperature to internal energy conversion factor */
- sp->temp_to_u_factor =
- phys_const->const_boltzmann_k /
- (p->mu_ionised * (hydro_gamma_minus_one)*phys_const->const_proton_mass);
-
- /* Fraction of energy in SNIa (?) */
- // Why is this one here? copied from EAGLE where it is always 1 ...
- sp->SNIa_energy_fraction = 1.0e0;
-
- /* Energy released by supernova */
- sp->total_energy_SNe =
- 1.0e51 / units_cgs_conversion_factor(us, UNIT_CONV_ENERGY);
-
- /* energy and temperature times h */
- sp->SNe_temperature =
- sp->total_energy_SNe / sp->temp_to_u_factor; // units_factor1 in EAGLE
-
- /* Find out timescale for feedback (used only for testing in const feedback
- * model) */
- sp->feedback_timescale =
- parser_get_opt_param_float(params, "Stars:feedback_timescale", 4e-5);
-
- /* Calculate number of supernovae per solar mass (used only for testing in
- * const feedback model) */
- const float ten_Myr_in_cgs = 3.154e14;
- const float SN_per_msun_factor = 0.01;
- sp->sn_per_msun =
- sp->feedback_timescale * units_cgs_conversion_factor(us, UNIT_CONV_TIME) /
- ten_Myr_in_cgs * SN_per_msun_factor; // timescale convert to cgs per 10
- // Myr (~3e14s). 0.01 solar masses
- // per supernova.
-
- /* Copy over solar mass (used only for testing in const feedback model) */
- sp->const_solar_mass = phys_const->const_solar_mass;
-}
-
-/**
- * @brief Print the global properties of the stars scheme.
- *
- * @param sp The #stars_props.
- */
-INLINE static void stars_props_print(const struct stars_props *sp) {
-
- /* Now stars */
- message("Stars kernel: %s with eta=%f (%.2f neighbours).", kernel_name,
- sp->eta_neighbours, sp->target_neighbours);
-
- message("Stars relative tolerance in h: %.5f (+/- %.4f neighbours).",
- sp->h_tolerance, sp->delta_neighbours);
-
- message(
- "Stars integration: Max change of volume: %.2f "
- "(max|dlog(h)/dt|=%f).",
- pow_dimension(expf(sp->log_max_h_change)), sp->log_max_h_change);
-
- if (sp->h_max != FLT_MAX)
- message("Maximal smoothing length allowed: %.4f", sp->h_max);
-
- message("Maximal iterations in ghost task set to %d",
- sp->max_smoothing_iterations);
-}
-
-#if defined(HAVE_HDF5)
-INLINE static void stars_props_print_snapshot(hid_t h_grpstars,
- const struct stars_props *sp) {
-
- io_write_attribute_s(h_grpstars, "Kernel function", kernel_name);
- io_write_attribute_f(h_grpstars, "Kernel target N_ngb",
- sp->target_neighbours);
- io_write_attribute_f(h_grpstars, "Kernel delta N_ngb", sp->delta_neighbours);
- io_write_attribute_f(h_grpstars, "Kernel eta", sp->eta_neighbours);
- io_write_attribute_f(h_grpstars, "Smoothing length tolerance",
- sp->h_tolerance);
- io_write_attribute_f(h_grpstars, "Maximal smoothing length", sp->h_max);
- io_write_attribute_f(h_grpstars, "Volume log(max(delta h))",
- sp->log_max_h_change);
- io_write_attribute_f(h_grpstars, "Volume max change time-step",
- pow_dimension(expf(sp->log_max_h_change)));
- io_write_attribute_i(h_grpstars, "Max ghost iterations",
- sp->max_smoothing_iterations);
-}
-#endif
-
-/**
- * @brief Write a #stars_props struct to the given FILE as a stream of bytes.
- *
- * @param p the struct
- * @param stream the file stream
- */
-INLINE static void stars_props_struct_dump(const struct stars_props *p,
- FILE *stream) {
- restart_write_blocks((void *)p, sizeof(struct stars_props), 1, stream,
- "starsprops", "stars props");
-}
-
-/**
- * @brief Restore a stars_props struct from the given FILE as a stream of
- * bytes.
- *
- * @param p the struct
- * @param stream the file stream
- */
-INLINE static void stars_props_struct_restore(const struct stars_props *p,
- FILE *stream) {
- restart_read_blocks((void *)p, sizeof(struct stars_props), 1, stream, NULL,
- "stars props");
-}
-
-#endif /* SWIFT_DEFAULT_STAR_IO_H */
diff --git a/src/stars/const/stars_part.h b/src/stars/const/stars_part.h
deleted file mode 100644
index 5927efff9b3abdb78dcdbf693d4044a53b1c73fa..0000000000000000000000000000000000000000
--- a/src/stars/const/stars_part.h
+++ /dev/null
@@ -1,189 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2016 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_DEFAULT_STAR_PART_H
-#define SWIFT_DEFAULT_STAR_PART_H
-
-/* Some standard headers. */
-#include <stdlib.h>
-
-/* Read chemistry */
-#include "chemistry_struct.h"
-
-/**
- * @brief Particle fields for the star particles.
- *
- * All quantities related to gravity are stored in the associate #gpart.
- */
-struct spart {
-
- /*! Particle ID. */
- long long id;
-
- /*! Pointer to corresponding gravity part. */
- struct gpart* gpart;
-
- /*! Particle position. */
- double x[3];
-
- /* Offset between current position and position at last tree rebuild. */
- float x_diff[3];
-
- /* Offset between current position and position at last tree rebuild. */
- float x_diff_sort[3];
-
- /*! Particle velocity. */
- float v[3];
-
- /*! Star mass */
- float mass;
-
- /*! Initial star mass */
- float mass_init;
-
- /* Particle cutoff radius. */
- float h;
-
- /*! Particle time bin */
- timebin_t time_bin;
-
- struct {
- /* Number of neighbours. */
- float wcount;
-
- /* Number of neighbours spatial derivative. */
- float wcount_dh;
-
- } density;
-
- struct {
- /* Change in smoothing length over time. */
- float h_dt;
-
- } feedback;
-
- struct {
- /* Mass of ejecta */
- float mass;
-
- /* Mass fractions of ejecta */
- struct chemistry_part_data chemistry_data;
-
- float ejecta_specific_thermal_energy;
-
- /* Number of type 1a SNe per unit mass */
- float num_SNIa;
-
- } to_distribute;
-
- /* kernel normalisation factor (equivalent to metalweight_norm in
- * eagle_enrich.c:811, TODO: IMPROVE COMMENT) */
- float omega_normalisation_inv;
-
- /* total mass of neighbouring gas particles */
- float ngb_mass;
-
- /*! Tracer structure */
- struct tracers_xpart_data tracers_data;
-
- /*! Chemistry structure */
- struct chemistry_part_data chemistry_data;
-
-#ifdef SWIFT_DEBUG_CHECKS
-
- /* Time of the last drift */
- integertime_t ti_drift;
-
- /* Time of the last kick */
- integertime_t ti_kick;
-
-#endif
-
-#ifdef DEBUG_INTERACTIONS_STARS
- /*! List of interacting particles in the density SELF and PAIR */
- long long ids_ngbs_density[MAX_NUM_OF_NEIGHBOURS_STARS];
-
- /*! Number of interactions in the density SELF and PAIR */
- int num_ngb_density;
-#endif
-
-} SWIFT_STRUCT_ALIGN;
-
-/**
- * @brief Contains all the constants and parameters of the stars scheme
- */
-struct stars_props {
-
- /*! Resolution parameter */
- float eta_neighbours;
-
- /*! Target weightd number of neighbours (for info only)*/
- float target_neighbours;
-
- /*! Smoothing length tolerance */
- float h_tolerance;
-
- /*! Tolerance on neighbour number (for info only)*/
- float delta_neighbours;
-
- /*! Maximal smoothing length */
- float h_max;
-
- /*! Maximal number of iterations to converge h */
- int max_smoothing_iterations;
-
- /*! Maximal change of h over one time-step */
- float log_max_h_change;
-
- /* Flag to switch between continuous and stochastic heating */
- int continuous_heating;
-
- /* Fraction of energy in SNIa (Note: always set to 1 in EAGLE, so may be not
- * necessary) */
- float SNIa_energy_fraction;
-
- /* Desired temperature increase due to supernovae */
- float SNe_deltaT_desired;
-
- /* Conversion factor from temperature to internal energy */
- float temp_to_u_factor;
-
- /* Energy released by one supernova */
- float total_energy_SNe;
-
- /* Temperature due to SNe (corresponding to units_factor1 in EAGLE) */
- float SNe_temperature;
-
- /* Timescale for feedback (used only for testing in const feedback model) */
- float feedback_timescale;
-
- /* Number of supernovae per solar mass (used only for testing in const
- * feedback model) */
- float sn_per_msun;
-
- /* Solar mass (used only for testing in const feedback model) */
- float const_solar_mass;
-
- /* Flag for testing energy injection */
- int const_feedback_energy_testing;
-
- // CHANGE THIS TO BE CONSISTENT WITH RAND MAX USED IN STAR FORMATION
- double inv_rand_max;
-};
-
-#endif /* SWIFT_DEFAULT_STAR_PART_H */