diff --git a/examples/main.c b/examples/main.c
index 2bfb7ea7600c119e067dc971fc494b74be0d9c44..bda2a0f57ad14e97262c5d77f3c9cc5de59b4559 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -912,10 +912,9 @@ int main(int argc, char *argv[]) {
bzero(&chemistry, sizeof(struct chemistry_global_data));
chemistry_init(params, &us, &prog_const, &chemistry);
if (myrank == 0) chemistry_print(&chemistry);
-
+
/* Initialise the cooling function properties */
- if (with_feedback)
- stars_evolve_init(params, &stars_properties);
+ if (with_feedback) stars_evolve_init(params, &stars_properties);
// ALEXEI: add print function!!!
/* Construct the engine policy */
diff --git a/src/cell.c b/src/cell.c
index f51c527dc18493b6b100e4ad3deca34df8c6143e..6ea2904ae919eda3be81a2842c61ca7bfb307ed0 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -3835,7 +3835,11 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
if (fabs(xp->v_full[0] * dt_drift) > e->s->dim[0] ||
fabs(xp->v_full[1] * dt_drift) > e->s->dim[1] ||
fabs(xp->v_full[2] * dt_drift) > e->s->dim[2]) {
- error("Particle drifts by more than a box length! id %llu xp->v_full %.5e %.5e %.5e p->v %.5e %.5e %.5e", p->id, xp->v_full[0], xp->v_full[1], xp->v_full[2], p->v[0], p->v[1], p->v[2]);
+ error(
+ "Particle drifts by more than a box length! id %llu xp->v_full "
+ "%.5e %.5e %.5e p->v %.5e %.5e %.5e",
+ p->id, xp->v_full[0], xp->v_full[1], xp->v_full[2], p->v[0],
+ p->v[1], p->v[2]);
}
#endif
@@ -3993,7 +3997,10 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
if (fabs(gp->v_full[0] * dt_drift) > e->s->dim[0] ||
fabs(gp->v_full[1] * dt_drift) > e->s->dim[1] ||
fabs(gp->v_full[2] * dt_drift) > e->s->dim[2]) {
- error("Particle drifts by more than a box length! gp->v_full %.5e %.5e %.5e", gp->v_full[0], gp->v_full[1], gp->v_full[2]);
+ error(
+ "Particle drifts by more than a box length! gp->v_full %.5e %.5e "
+ "%.5e",
+ gp->v_full[0], gp->v_full[1], gp->v_full[2]);
}
#endif
diff --git a/src/engine.c b/src/engine.c
index 43a0e96972333d4bfb2f6f3130f9dd6bb387cc6a..9f8f564f5bdd0b671651cb2b5ddc4c70e2cb68dc 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2882,7 +2882,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
if (!flag_entropy_ICs) {
if (e->nodeID == 0) message("Converting internal energy variable.");
-
+
space_convert_quantities(e->s, e->verbose);
/* Correct what we did (e.g. in PE-SPH, need to recompute rho_bar) */
@@ -2913,7 +2913,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
/* Prepare all the tasks again for a new round */
engine_marktasks(e);
-
+
/* No drift this time */
engine_skip_drift(e);
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 1a513e5dada59c1d776d85b71e77fccd4056f101..6793698f4484ce74cdb3c86a55069cc172cb0e61 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -360,8 +360,7 @@ __attribute__((always_inline)) INLINE static void
hydro_set_physical_internal_energy(struct part *restrict p,
const struct cosmology *restrict cosmo,
float u) {
- p->entropy =
- gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, u);
+ p->entropy = gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, u);
}
/**
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index ba186f36ab57c47c77193c189d1a242e87dd0d9c..80a7787f111c83debffb8f0e33e3c902fc9f39c3 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -527,16 +527,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Eventually got the acceleration */
const float acc = visc_term + sph_term;
- //if (pi->id == 142820 || pi->id == 145267) message("1 id %llu acc %.5e visc_term %.5e sph_term %.5e", pi->id, acc, visc_term, sph_term);
- //if (pj->id == 142820 || pj->id == 145267) message("2 id %llu acc %.5e visc_term %.5e sph_term %.5e", pj->id, acc, visc_term, sph_term);
+ // if (pi->id == 142820 || pi->id == 145267) message("1 id %llu acc %.5e
+ // visc_term %.5e sph_term %.5e", pi->id, acc, visc_term, sph_term); if (pj->id
+ // == 142820 || pj->id == 145267) message("2 id %llu acc %.5e visc_term %.5e
+ // sph_term %.5e", pj->id, acc, visc_term, sph_term);
/* Use the force Luke ! */
- //if (pi->id == 145267) message("pi 1 id %llu acted on by id %llu a_old %.5e %.5e %.5e mj %.5e acc %.5e dx %.5e %.5e %.5e a_new %.5e %.5e %.5e", pi->id, pj->id, pi->a_hydro[0], pi->a_hydro[1], pi->a_hydro[2], mj, acc, dx[0], dx[1], dx[2], pi->a_hydro[0] + mj * acc * dx[0], pi->a_hydro[1] + mj * acc * dx[1], pi->a_hydro[2] + mj * acc * dx[2]);
+ // if (pi->id == 145267) message("pi 1 id %llu acted on by id %llu a_old %.5e
+ // %.5e %.5e mj %.5e acc %.5e dx %.5e %.5e %.5e a_new %.5e %.5e %.5e", pi->id,
+ // pj->id, pi->a_hydro[0], pi->a_hydro[1], pi->a_hydro[2], mj, acc, dx[0],
+ // dx[1], dx[2], pi->a_hydro[0] + mj * acc * dx[0], pi->a_hydro[1] + mj * acc *
+ // dx[1], pi->a_hydro[2] + mj * acc * dx[2]);
pi->a_hydro[0] -= mj * acc * dx[0];
pi->a_hydro[1] -= mj * acc * dx[1];
pi->a_hydro[2] -= mj * acc * dx[2];
- //if (pj->id == 145267) message("pj id %llu acted on by %llu a_old %.5e %.5e %.5e mi %.5e acc %.5e dx %.5e %.5e %.5e a_new %.5e %.5e %.5e", pj->id, pi->id, pj->a_hydro[0], pj->a_hydro[1], pj->a_hydro[2], mi, acc, dx[0], dx[1], dx[2], pj->a_hydro[0] + mi * acc * dx[0], pj->a_hydro[1] + mi * acc * dx[1], pj->a_hydro[2] + mi * acc * dx[2]);
+ // if (pj->id == 145267) message("pj id %llu acted on by %llu a_old %.5e %.5e
+ // %.5e mi %.5e acc %.5e dx %.5e %.5e %.5e a_new %.5e %.5e %.5e", pj->id,
+ // pi->id, pj->a_hydro[0], pj->a_hydro[1], pj->a_hydro[2], mi, acc, dx[0],
+ // dx[1], dx[2], pj->a_hydro[0] + mi * acc * dx[0], pj->a_hydro[1] + mi * acc *
+ // dx[1], pj->a_hydro[2] + mi * acc * dx[2]);
pj->a_hydro[0] += mi * acc * dx[0];
pj->a_hydro[1] += mi * acc * dx[1];
pj->a_hydro[2] += mi * acc * dx[2];
@@ -660,10 +670,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Eventually got the acceleration */
const float acc = visc_term + sph_term;
- //if (pi->id == 142820 || pi->id == 145267) message("3 id %llu acc %.5e visc_term %.5e sph_term %.5e", pi->id, acc, visc_term, sph_term);
+ // if (pi->id == 142820 || pi->id == 145267) message("3 id %llu acc %.5e
+ // visc_term %.5e sph_term %.5e", pi->id, acc, visc_term, sph_term);
/* Use the force Luke ! */
- //if (pi->id == 145267) message("pi 2 id %llu acted on by %llu a_old %.5e %.5e %.5e mj %.5e acc %.5e dx %.5e %.5e %.5e a_new %.5e %.5e %.5e", pi->id, pj->id, pi->a_hydro[0], pi->a_hydro[1], pi->a_hydro[2], mj, acc, dx[0], dx[1], dx[2], pi->a_hydro[0] + mj * acc * dx[0], pi->a_hydro[1] + mj * acc * dx[1], pi->a_hydro[2] + mj * acc * dx[2]);
+ // if (pi->id == 145267) message("pi 2 id %llu acted on by %llu a_old %.5e
+ // %.5e %.5e mj %.5e acc %.5e dx %.5e %.5e %.5e a_new %.5e %.5e %.5e", pi->id,
+ // pj->id, pi->a_hydro[0], pi->a_hydro[1], pi->a_hydro[2], mj, acc, dx[0],
+ // dx[1], dx[2], pi->a_hydro[0] + mj * acc * dx[0], pi->a_hydro[1] + mj * acc *
+ // dx[1], pi->a_hydro[2] + mj * acc * dx[2]);
pi->a_hydro[0] -= mj * acc * dx[0];
pi->a_hydro[1] -= mj * acc * dx[1];
pi->a_hydro[2] -= mj * acc * dx[2];
@@ -1077,16 +1092,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_limiter(
/* Wake up the neighbour? */
if (pi->force.v_sig > const_limiter_max_v_sig_ratio * pj->force.v_sig) {
- // ALEXEI seems to crash with this option when running with debug checks
+ // ALEXEI seems to crash with this option when running with debug checks
// on comparing ti_kick to ti_start in kick_part. Use code below instead
// (commented MATTHIEU)
- //pj->wakeup = time_bin_awake;
+ // pj->wakeup = time_bin_awake;
// MATTHIEU
- if (pj->wakeup == time_bin_not_awake)
- pj->wakeup = time_bin_awake;
- else if (pj->wakeup > 0)
- pj->wakeup = -pj->wakeup;
+ if (pj->wakeup == time_bin_not_awake)
+ pj->wakeup = time_bin_awake;
+ else if (pj->wakeup > 0)
+ pj->wakeup = -pj->wakeup;
}
}
diff --git a/src/runner.c b/src/runner.c
index a2678bc89cebff1a24adb5a71c71edcc5d3f3219..a3dd975210a9f85a81ebf6a4f04afeb52f38c9f5 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -179,7 +179,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
/* While there are particles that need to be updated... */
for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
num_reruns++) {
-
+
/* Reset the redo-count. */
redo = 0;
@@ -189,14 +189,14 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
/* Get a direct pointer on the part. */
struct spart *sp = &sparts[sid[i]];
- /* get particle timestep */
+ /* get particle timestep */
double dt;
if (with_cosmology) {
const integertime_t ti_step = get_integer_timestep(sp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(e->ti_current - 1, sp->time_bin);
dt = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
- ti_begin + ti_step);
+ ti_begin + ti_step);
} else {
dt = get_timestep(sp->time_bin, e->time_base);
}
@@ -351,11 +351,14 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
/* We now have a particle whose smoothing length has converged */
stars_reset_feedback(sp);
- // Get current time
- float current_time_begin = get_integer_time_begin(e->ti_current - 1, sp->time_bin) * e->time_base;
+ // Get current time
+ float current_time_begin =
+ get_integer_time_begin(e->ti_current - 1, sp->time_bin) *
+ e->time_base;
/* Compute the stellar evolution */
- stars_evolve_spart(sp, e->stars_properties, cosmo, us, current_time_begin, dt);
+ stars_evolve_spart(sp, e->stars_properties, cosmo, us,
+ current_time_begin, dt);
}
/* We now need to treat the particles whose smoothing length had not
diff --git a/src/runner_doiact_stars.h b/src/runner_doiact_stars.h
index 66b947af6b6e7ea156983ce2c6ac014f2a926099..2a45fce24eb085bf4b9681409c66c88653efc8c8 100644
--- a/src/runner_doiact_stars.h
+++ b/src/runner_doiact_stars.h
@@ -128,7 +128,8 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
#endif
if (r2 > 0.f && r2 < hig2) {
- IACT_STARS(r2, dx, hi, hj, si, pj, cosmo, stars_properties, xpj, e->ti_current);
+ IACT_STARS(r2, dx, hi, hj, si, pj, cosmo, stars_properties, xpj,
+ e->ti_current);
}
} /* loop over the parts in ci. */
} /* loop over the sparts in ci. */
@@ -207,7 +208,9 @@ void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
error("Particle pj not drifted to current time");
#endif
- if (r2 < hig2) IACT_STARS(r2, dx, hi, hj, si, pj, cosmo, stars_properties, xpj, e->ti_current);
+ if (r2 < hig2)
+ IACT_STARS(r2, dx, hi, hj, si, pj, cosmo, stars_properties, xpj,
+ e->ti_current);
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -297,7 +300,8 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
#endif
/* Hit or miss? */
if (r2 < hig2) {
- IACT_STARS(r2, dx, hi, pj->h, spi, pj, cosmo, stars_properties, xpj, e->ti_current);
+ IACT_STARS(r2, dx, hi, pj->h, spi, pj, cosmo, stars_properties, xpj,
+ e->ti_current);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -368,7 +372,8 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
/* Hit or miss? */
if (r2 > 0.f && r2 < hig2) {
- IACT_STARS(r2, dx, hi, hj, spi, pj, cosmo, stars_properties, xpj, e->ti_current);
+ IACT_STARS(r2, dx, hi, hj, spi, pj, cosmo, stars_properties, xpj,
+ e->ti_current);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
diff --git a/src/stars/Default/stars_iact.h b/src/stars/Default/stars_iact.h
index 2b5154c0f32c38ea1d64804a25a5b11eeaf0db82..8c83807f1b97b53e62025007d3d6f05ac3dd3755 100644
--- a/src/stars/Default/stars_iact.h
+++ b/src/stars/Default/stars_iact.h
@@ -33,12 +33,12 @@
* @param H Current Hubble parameter.
*/
__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, float a,
- float H, const struct cosmology *restrict cosmo,
- const struct stars_props *restrict stars_properties,
- struct xpart *restrict xp) {
+runner_iact_nonsym_stars_density(
+ float r2, const float *dx, float hi, float hj, struct spart *restrict si,
+ const struct part *restrict pj, float a, float H,
+ const struct cosmology *restrict cosmo,
+ const struct stars_props *restrict stars_properties,
+ struct xpart *restrict xp) {
float wi, wi_dx;
@@ -78,12 +78,12 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
* @param H Current Hubble parameter.
*/
__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, float a, float H,
- const struct cosmology *restrict cosmo,
- const struct stars_props *restrict stars_properties,
- struct xpart *restrict xp) {
+runner_iact_nonsym_stars_feedback(
+ float r2, const float *dx, float hi, float hj, struct spart *restrict si,
+ struct part *restrict pj, float a, float H,
+ const struct cosmology *restrict cosmo,
+ const struct stars_props *restrict stars_properties,
+ struct xpart *restrict xp) {
const float mj = hydro_get_mass(pj);
const float rhoj = pj->rho;
diff --git a/src/stars/Default/stars_io.h b/src/stars/Default/stars_io.h
index 2cff7b604f7a949c83b0a09e6a5c917776dbd2ec..ce54191e39901a1056abdefecb069563e60fd0cb 100644
--- a/src/stars/Default/stars_io.h
+++ b/src/stars/Default/stars_io.h
@@ -109,7 +109,7 @@ INLINE static void stars_props_init(struct stars_props *sp,
const struct unit_system *us,
struct swift_params *params,
const struct hydro_props *p,
- const struct cosmology *cosmo) {
+ const struct cosmology *cosmo) {
/* Kernel properties */
sp->eta_neighbours = parser_get_opt_param_float(
diff --git a/src/stars/EAGLE/imf.h b/src/stars/EAGLE/imf.h
index fef12d03ca18f31bdb287910df88229243a89cb7..7fe8c69c70d90443419591a15b8c6e88da9d6057 100644
--- a/src/stars/EAGLE/imf.h
+++ b/src/stars/EAGLE/imf.h
@@ -27,38 +27,46 @@ static const int N_imf_mass_bins = 200;
static const int gamma_SNIa = 2;
// do we need doubles in signature?
-inline static void determine_imf_bins(double log_min_dying_mass, double log_max_dying_mass,
- int *ilow, int *ihigh,
- const struct stars_props *restrict star_properties) {
+inline static void determine_imf_bins(
+ double log_min_dying_mass, double log_max_dying_mass, int *ilow, int *ihigh,
+ const struct stars_props *restrict star_properties) {
int i1, i2;
if (log_min_dying_mass < star_properties->imf_mass_bin_log10[0])
log_min_dying_mass = star_properties->imf_mass_bin_log10[0];
- if (log_min_dying_mass > star_properties->imf_mass_bin_log10[N_imf_mass_bins - 1])
- log_min_dying_mass = star_properties->imf_mass_bin_log10[N_imf_mass_bins - 1];
+ if (log_min_dying_mass >
+ star_properties->imf_mass_bin_log10[N_imf_mass_bins - 1])
+ log_min_dying_mass =
+ star_properties->imf_mass_bin_log10[N_imf_mass_bins - 1];
if (log_max_dying_mass < star_properties->imf_mass_bin_log10[0])
log_max_dying_mass = star_properties->imf_mass_bin_log10[0];
- if (log_max_dying_mass > star_properties->imf_mass_bin_log10[N_imf_mass_bins - 1])
- log_max_dying_mass = star_properties->imf_mass_bin_log10[N_imf_mass_bins - 1];
+ if (log_max_dying_mass >
+ star_properties->imf_mass_bin_log10[N_imf_mass_bins - 1])
+ log_max_dying_mass =
+ star_properties->imf_mass_bin_log10[N_imf_mass_bins - 1];
- for (i1 = 0;
- i1 < N_imf_mass_bins - 2 && star_properties->imf_mass_bin_log10[i1 + 1] < log_min_dying_mass;
- i1++);
-
- for (i2 = 1;
- i2 < N_imf_mass_bins - 1 && star_properties->imf_mass_bin_log10[i2] < log_max_dying_mass; i2++);
+ for (i1 = 0; i1 < N_imf_mass_bins - 2 &&
+ star_properties->imf_mass_bin_log10[i1 + 1] < log_min_dying_mass;
+ i1++)
+ ;
+ for (i2 = 1; i2 < N_imf_mass_bins - 1 &&
+ star_properties->imf_mass_bin_log10[i2] < log_max_dying_mass;
+ i2++)
+ ;
*ilow = i1;
*ihigh = i2;
}
-// Change this function to not pass in stellar_yields as that is in star_properties.
-inline static float integrate_imf(float log_min_mass, float log_max_mass, float m2, int mode, float *stellar_yields,
- const struct stars_props *restrict star_properties){
+// Change this function to not pass in stellar_yields as that is in
+// star_properties.
+inline static float integrate_imf(
+ float log_min_mass, float log_max_mass, float m2, int mode,
+ float *stellar_yields, const struct stars_props *restrict star_properties) {
double result, u, f;
@@ -66,32 +74,39 @@ inline static float integrate_imf(float log_min_mass, float log_max_mass, float
float dlm, dm;
- dlm = star_properties->imf_mass_bin_log10[1] - star_properties->imf_mass_bin_log10[0]; /* dlog(m) */
+ dlm = star_properties->imf_mass_bin_log10[1] -
+ star_properties->imf_mass_bin_log10[0]; /* dlog(m) */
- determine_imf_bins(log_min_mass, log_max_mass, &ilow, &ihigh, star_properties);
+ determine_imf_bins(log_min_mass, log_max_mass, &ilow, &ihigh,
+ star_properties);
float integrand[N_imf_mass_bins];
if (mode == 0)
for (index = ilow; index < ihigh + 1; index++)
integrand[index] =
- star_properties->imf_by_number[index] * star_properties->imf_mass_bin[index]; /* integrate by number */
+ star_properties->imf_by_number[index] *
+ star_properties->imf_mass_bin[index]; /* integrate by number */
else if (mode == 1)
for (index = ilow; index < ihigh + 1; index++)
- integrand[index] = star_properties->imf_by_number[index] * star_properties->imf_mass_bin[index] *
- star_properties->imf_mass_bin[index]; /* integrate by mass */
+ integrand[index] =
+ star_properties->imf_by_number[index] *
+ star_properties->imf_mass_bin[index] *
+ star_properties->imf_mass_bin[index]; /* integrate by mass */
else if (mode == 2)
for (index = ilow; index < ihigh + 1; index++)
integrand[index] =
stellar_yields[index] * star_properties->imf_by_number[index] *
- star_properties->imf_mass_bin[index]; /* integrate number * yield weighted */
+ star_properties
+ ->imf_mass_bin[index]; /* integrate number * yield weighted */
else if (mode == 3)
for (index = ilow; index < ihigh + 1; index++) {
u = m2 / star_properties->imf_mass_bin[index];
f = pow(2.0, gamma_SNIa + 1) * (gamma_SNIa + 1) * pow(u, gamma_SNIa);
integrand[index] =
f * star_properties->imf_by_number[index] /
- star_properties->imf_mass_bin[index]; /* integrate number * f(u) / M ... type Ia SN */
+ star_properties->imf_mass_bin[index]; /* integrate number * f(u) / M
+ ... type Ia SN */
}
else {
error("invalid mode in integrate_imf = %d\n", mode);
@@ -127,24 +142,34 @@ inline static float integrate_imf(float log_min_mass, float log_max_mass, float
function of log_10(mass) */
return result;
-
}
-inline static void init_imf(struct stars_props *restrict star_properties){
+inline static void init_imf(struct stars_props *restrict star_properties) {
// ALEXEI: use better names for solar_mass, log_solar_mass
float norm = 0, solar_mass, log_solar_mass;
- const float dlm = (log_imf_max_solar_mass - log_imf_min_solar_mass) / (double)(N_imf_mass_bins - 1);
+ const float dlm = (log_imf_max_solar_mass - log_imf_min_solar_mass) /
+ (double)(N_imf_mass_bins - 1);
- if (posix_memalign((void **)&star_properties->imf_by_number, SWIFT_STRUCT_ALIGNMENT, N_imf_mass_bins*sizeof(float)) != 0)
+ if (posix_memalign((void **)&star_properties->imf_by_number,
+ SWIFT_STRUCT_ALIGNMENT,
+ N_imf_mass_bins * sizeof(float)) != 0)
error("Failed to allocate IMF bins table");
- if (posix_memalign((void **)&star_properties->imf_mass_bin, SWIFT_STRUCT_ALIGNMENT, N_imf_mass_bins*sizeof(float)) != 0)
+ if (posix_memalign((void **)&star_properties->imf_mass_bin,
+ SWIFT_STRUCT_ALIGNMENT,
+ N_imf_mass_bins * sizeof(float)) != 0)
error("Failed to allocate IMF bins table");
- if (posix_memalign((void **)&star_properties->imf_mass_bin_log10, SWIFT_STRUCT_ALIGNMENT, N_imf_mass_bins*sizeof(float)) != 0)
+ if (posix_memalign((void **)&star_properties->imf_mass_bin_log10,
+ SWIFT_STRUCT_ALIGNMENT,
+ N_imf_mass_bins * sizeof(float)) != 0)
error("Failed to allocate IMF bins table");
- if (posix_memalign((void **)&star_properties->imf_by_number1, SWIFT_STRUCT_ALIGNMENT, N_imf_mass_bins*sizeof(float)) != 0)
+ if (posix_memalign((void **)&star_properties->imf_by_number1,
+ SWIFT_STRUCT_ALIGNMENT,
+ N_imf_mass_bins * sizeof(float)) != 0)
error("Failed to allocate IMF bins table");
- if (posix_memalign((void **)&star_properties->stellar_yield, SWIFT_STRUCT_ALIGNMENT, N_imf_mass_bins*sizeof(float)) != 0)
+ if (posix_memalign((void **)&star_properties->stellar_yield,
+ SWIFT_STRUCT_ALIGNMENT,
+ N_imf_mass_bins * sizeof(float)) != 0)
error("Failed to allocate IMF bins table");
float dummy_stellar_fields;
@@ -160,7 +185,8 @@ inline static void init_imf(struct stars_props *restrict star_properties){
solar_mass = exp(M_LN10 * log_solar_mass);
// can these pows be replaced with some trick?
- star_properties->imf_by_number[i] = pow(solar_mass, -star_properties->IMF_Exponent);
+ star_properties->imf_by_number[i] =
+ pow(solar_mass, -star_properties->IMF_Exponent);
star_properties->imf_mass_bin[i] = solar_mass;
star_properties->imf_mass_bin_log10[i] = log_solar_mass;
@@ -176,23 +202,25 @@ inline static void init_imf(struct stars_props *restrict star_properties){
if (solar_mass > 1.0)
star_properties->imf_by_number[i] = 0.237912 * pow(solar_mass, -2.3);
else
- star_properties->imf_by_number[i] = 0.852464 *
- exp((log10(solar_mass) - log10(0.079)) * (log10(solar_mass) - log10(0.079)) /
- (-2.0 * pow(0.69, 2))) /
- solar_mass;
+ star_properties->imf_by_number[i] =
+ 0.852464 *
+ exp((log10(solar_mass) - log10(0.079)) *
+ (log10(solar_mass) - log10(0.079)) / (-2.0 * pow(0.69, 2))) /
+ solar_mass;
star_properties->imf_mass_bin[i] = solar_mass;
star_properties->imf_mass_bin_log10[i] = log_solar_mass;
}
} else {
- error("Invalid IMF model %s. Valid models are: PowerLaw and Chabrier\n", star_properties->IMF_Model);
+ error("Invalid IMF model %s. Valid models are: PowerLaw and Chabrier\n",
+ star_properties->IMF_Model);
}
- norm = integrate_imf(log_imf_min_solar_mass, log_imf_max_solar_mass,
- 0.0, 1, &dummy_stellar_fields, star_properties);
-
- for (int i = 0; i < N_imf_mass_bins; i++) star_properties->imf_by_number[i] /= norm;
+ norm = integrate_imf(log_imf_min_solar_mass, log_imf_max_solar_mass, 0.0, 1,
+ &dummy_stellar_fields, star_properties);
+ for (int i = 0; i < N_imf_mass_bins; i++)
+ star_properties->imf_by_number[i] /= norm;
}
#endif
diff --git a/src/stars/EAGLE/stars.h b/src/stars/EAGLE/stars.h
index c0cfb89020b652829471e8eb4c3af70f2a3b2ddc..adc7f69c76c21adeb2801e796feb3bc0e6197387 100644
--- a/src/stars/EAGLE/stars.h
+++ b/src/stars/EAGLE/stars.h
@@ -20,13 +20,13 @@
#define SWIFT_EAGLE_STARS_H
#include <float.h>
-#include "minmax.h"
#include "imf.h"
+#include "minmax.h"
#include "yield_tables.h"
-static const float log10_SNII_min_mass_msun = 0.77815125f; // log10(6).
-static const float log10_SNII_max_mass_msun = 2.f; // log10(100).
-static const float log10_SNIa_max_mass_msun = 0.90308999f; // log10(8).
+static const float log10_SNII_min_mass_msun = 0.77815125f; // log10(6).
+static const float log10_SNII_max_mass_msun = 2.f; // log10(100).
+static const float log10_SNIa_max_mass_msun = 0.90308999f; // log10(8).
static const float imf_max_mass_msun = 100.f;
static const float imf_min_mass_msun = 0.1;
@@ -79,13 +79,13 @@ __attribute__((always_inline)) INLINE static void stars_first_init_spart(
sp->birth_time = -1.f;
// ALEXEI: initialise chemistry data for feedback tests
- for (int i = 0; i < chemistry_element_count; i++) sp->chemistry_data.metal_mass_fraction[i] = 0.f;
+ for (int i = 0; i < chemistry_element_count; i++)
+ sp->chemistry_data.metal_mass_fraction[i] = 0.f;
sp->chemistry_data.metal_mass_fraction[0] = 0.742;
sp->chemistry_data.metal_mass_fraction[1] = 0.238;
sp->chemistry_data.metal_mass_fraction_total = 0.02;
-
stars_init_spart(sp);
}
@@ -195,7 +195,7 @@ __attribute__((always_inline)) INLINE static void stars_reset_acceleration(
// -------------------- Work in progress ------------------------------
// This really needs to be somewhere else
-inline static double interpol_1d(double *table, int i, float dx) {
+inline static double interpol_1d(double* table, int i, float dx) {
double result;
result = (1 - dx) * table[i] + dx * table[i + 1];
@@ -203,7 +203,8 @@ inline static double interpol_1d(double *table, int i, float dx) {
return result;
}
-inline static double interpol_2d(double **table, int i, int j, float dx, float dy) {
+inline static double interpol_2d(double** table, int i, int j, float dx,
+ float dy) {
double result;
result = (1 - dx) * (1 - dy) * table[i][j] + (1 - dx) * dy * table[i][j + 1] +
@@ -212,10 +213,9 @@ inline static double interpol_2d(double **table, int i, int j, float dx, float d
return result;
}
-
-
inline static float dying_mass_msun(float age_Gyr, float metallicity,
- const struct stars_props *restrict star_properties) {
+ const struct stars_props* restrict
+ star_properties) {
// check out units for all these quantities and name them accordingly
float mass = 0, d_metal, d_time1 = 0, d_time2 = 0, log_age_yr, mass1, mass2;
@@ -223,22 +223,22 @@ inline static float dying_mass_msun(float age_Gyr, float metallicity,
int metal_index, i, index_time1 = -1, index_time2 = -1;
// What do we do with the constants here?
- switch(star_properties->stellar_lifetime_flag) {
+ switch (star_properties->stellar_lifetime_flag) {
case 0:
/* padovani & matteucci 1993 */
if (age_Gyr > 0.039765318659064693)
- mass = exp(M_LN10 *
- (7.764 - (1.79 -
- (1.338 - 0.1116 * (9 + log10(age_Gyr))) *
- (1.338 - 0.1116 * (9 + log10(age_Gyr)))) /
- 0.2232));
+ mass =
+ exp(M_LN10 *
+ (7.764 - (1.79 - (1.338 - 0.1116 * (9 + log10(age_Gyr))) *
+ (1.338 - 0.1116 * (9 + log10(age_Gyr)))) /
+ 0.2232));
else if (age_Gyr > 0.003)
mass = pow((age_Gyr - 0.003) / 1.2, -1 / 1.85);
else
mass = imf_max_mass_msun;
break;
-
+
case 1:
/* maeder & meynet 1989 */
@@ -272,14 +272,20 @@ inline static float dying_mass_msun(float age_Gyr, float metallicity,
if (metallicity <= star_properties->lifetimes.metallicity[0]) {
metal_index = 0;
d_metal = 0.0;
- } else if (metallicity >= star_properties->lifetimes.metallicity[star_properties->lifetimes.n_z - 1]) {
+ } else if (metallicity >=
+ star_properties->lifetimes
+ .metallicity[star_properties->lifetimes.n_z - 1]) {
metal_index = star_properties->lifetimes.n_z - 2;
d_metal = 1.0;
} else {
- for (metal_index = 0; metal_index < star_properties->lifetimes.n_z - 1; metal_index++)
- if (star_properties->lifetimes.metallicity[metal_index + 1] > metallicity) break;
-
- d_metal = (metallicity - star_properties->lifetimes.metallicity[metal_index]) /
+ for (metal_index = 0; metal_index < star_properties->lifetimes.n_z - 1;
+ metal_index++)
+ if (star_properties->lifetimes.metallicity[metal_index + 1] >
+ metallicity)
+ break;
+
+ d_metal = (metallicity -
+ star_properties->lifetimes.metallicity[metal_index]) /
(star_properties->lifetimes.metallicity[metal_index + 1] -
star_properties->lifetimes.metallicity[metal_index]);
}
@@ -287,16 +293,22 @@ inline static float dying_mass_msun(float age_Gyr, float metallicity,
if (log_age_yr >= star_properties->lifetimes.dyingtime[metal_index][0]) {
index_time1 = 0;
d_time1 = 0.0;
- } else if (log_age_yr <= star_properties->lifetimes.dyingtime[metal_index][star_properties->lifetimes.n_mass - 1]) {
+ } else if (log_age_yr <=
+ star_properties->lifetimes
+ .dyingtime[metal_index]
+ [star_properties->lifetimes.n_mass - 1]) {
index_time1 = star_properties->lifetimes.n_mass - 2;
d_time1 = 1.0;
}
- if (log_age_yr >= star_properties->lifetimes.dyingtime[metal_index + 1][0]) {
+ if (log_age_yr >=
+ star_properties->lifetimes.dyingtime[metal_index + 1][0]) {
index_time2 = 0;
d_time2 = 0.0;
} else if (log_age_yr <=
- star_properties->lifetimes.dyingtime[metal_index + 1][star_properties->lifetimes.n_mass - 1]) {
+ star_properties->lifetimes
+ .dyingtime[metal_index + 1]
+ [star_properties->lifetimes.n_mass - 1]) {
index_time2 = star_properties->lifetimes.n_mass - 2;
d_time2 = 1.0;
}
@@ -304,29 +316,41 @@ inline static float dying_mass_msun(float age_Gyr, float metallicity,
i = star_properties->lifetimes.n_mass;
while (i >= 0 && (index_time1 == -1 || index_time2 == -1)) {
i--;
- if (star_properties->lifetimes.dyingtime[metal_index][i] >= log_age_yr && index_time1 == -1) {
+ if (star_properties->lifetimes.dyingtime[metal_index][i] >=
+ log_age_yr &&
+ index_time1 == -1) {
index_time1 = i;
- d_time1 = (log_age_yr - star_properties->lifetimes.dyingtime[metal_index][index_time1]) /
- (star_properties->lifetimes.dyingtime[metal_index][index_time1 + 1] -
- star_properties->lifetimes.dyingtime[metal_index][index_time1]);
+ d_time1 =
+ (log_age_yr -
+ star_properties->lifetimes.dyingtime[metal_index][index_time1]) /
+ (star_properties->lifetimes
+ .dyingtime[metal_index][index_time1 + 1] -
+ star_properties->lifetimes.dyingtime[metal_index][index_time1]);
}
- if (star_properties->lifetimes.dyingtime[metal_index + 1][i] >= log_age_yr && index_time2 == -1) {
+ if (star_properties->lifetimes.dyingtime[metal_index + 1][i] >=
+ log_age_yr &&
+ index_time2 == -1) {
index_time2 = i;
- d_time2 = (log_age_yr - star_properties->lifetimes.dyingtime[metal_index + 1][index_time2]) /
- (star_properties->lifetimes.dyingtime[metal_index + 1][index_time2 + 1] -
- star_properties->lifetimes.dyingtime[metal_index + 1][index_time2]);
+ d_time2 =
+ (log_age_yr - star_properties->lifetimes
+ .dyingtime[metal_index + 1][index_time2]) /
+ (star_properties->lifetimes
+ .dyingtime[metal_index + 1][index_time2 + 1] -
+ star_properties->lifetimes
+ .dyingtime[metal_index + 1][index_time2]);
}
}
- mass1 = interpol_1d(star_properties->lifetimes.mass, index_time1, d_time1);
- mass2 = interpol_1d(star_properties->lifetimes.mass, index_time2, d_time2);
+ mass1 =
+ interpol_1d(star_properties->lifetimes.mass, index_time1, d_time1);
+ mass2 =
+ interpol_1d(star_properties->lifetimes.mass, index_time2, d_time2);
mass = (1.0 - d_metal) * mass1 + d_metal * mass2;
break;
default:
error("stellar lifetimes not defined\n");
-
}
if (mass > imf_max_mass_msun) mass = imf_max_mass_msun;
@@ -334,7 +358,8 @@ inline static float dying_mass_msun(float age_Gyr, float metallicity,
}
inline static float lifetime_in_Gyr(float mass, float metallicity,
- const struct stars_props *restrict star_properties) {
+ const struct stars_props* restrict
+ star_properties) {
double time = 0, d_mass, d_metal;
@@ -377,35 +402,46 @@ inline static float lifetime_in_Gyr(float mass, float metallicity,
if (mass <= star_properties->lifetimes.mass[0]) {
mass_index = 0;
d_mass = 0.0;
- } else if (mass >= star_properties->lifetimes.mass[star_properties->lifetimes.n_mass - 1]) {
+ } else if (mass >= star_properties->lifetimes
+ .mass[star_properties->lifetimes.n_mass - 1]) {
mass_index = star_properties->lifetimes.n_mass - 2;
d_mass = 1.0;
} else {
- for (mass_index = 0; mass_index < star_properties->lifetimes.n_mass - 1; mass_index++)
+ for (mass_index = 0; mass_index < star_properties->lifetimes.n_mass - 1;
+ mass_index++)
if (star_properties->lifetimes.mass[mass_index + 1] > mass) break;
d_mass = (mass - star_properties->lifetimes.mass[mass_index]) /
- (star_properties->lifetimes.mass[mass_index + 1] - star_properties->lifetimes.mass[mass_index]);
+ (star_properties->lifetimes.mass[mass_index + 1] -
+ star_properties->lifetimes.mass[mass_index]);
}
if (metallicity <= star_properties->lifetimes.metallicity[0]) {
metal_index = 0;
d_metal = 0.0;
- } else if (metallicity >= star_properties->lifetimes.metallicity[star_properties->lifetimes.n_z - 1]) {
+ } else if (metallicity >=
+ star_properties->lifetimes
+ .metallicity[star_properties->lifetimes.n_z - 1]) {
metal_index = star_properties->lifetimes.n_z - 2;
d_metal = 1.0;
} else {
- for (metal_index = 0; metal_index < star_properties->lifetimes.n_z - 1; metal_index++)
- if (star_properties->lifetimes.metallicity[metal_index + 1] > metallicity) break;
-
- d_metal = (metallicity - star_properties->lifetimes.metallicity[metal_index]) /
+ for (metal_index = 0; metal_index < star_properties->lifetimes.n_z - 1;
+ metal_index++)
+ if (star_properties->lifetimes.metallicity[metal_index + 1] >
+ metallicity)
+ break;
+
+ d_metal = (metallicity -
+ star_properties->lifetimes.metallicity[metal_index]) /
(star_properties->lifetimes.metallicity[metal_index + 1] -
star_properties->lifetimes.metallicity[metal_index]);
}
/* time in Gyr */
- time = exp(M_LN10 * interpol_2d(star_properties->lifetimes.dyingtime, metal_index, mass_index,
- d_metal, d_mass)) / 1.0e9;
+ time =
+ exp(M_LN10 * interpol_2d(star_properties->lifetimes.dyingtime,
+ metal_index, mass_index, d_metal, d_mass)) /
+ 1.0e9;
break;
@@ -416,26 +452,31 @@ inline static float lifetime_in_Gyr(float mass, float metallicity,
return time;
}
-inline static void determine_bin_yield_AGB(int *iz_low, int *iz_high, float *dz, float log_metallicity,
- const struct stars_props *restrict star_properties){
+inline static void determine_bin_yield_AGB(
+ int* iz_low, int* iz_high, float* dz, float log_metallicity,
+ const struct stars_props* restrict star_properties) {
if (log_metallicity > log_min_metallicity) {
int j;
for (j = 0; j < star_properties->AGB_n_z - 1 &&
- log_metallicity > star_properties->yield_AGB.metallicity[j + 1];
- j++);
+ log_metallicity > star_properties->yield_AGB.metallicity[j + 1];
+ j++)
+ ;
*iz_low = j;
*iz_high = *iz_low + 1;
- if (*iz_high >= star_properties->AGB_n_z) *iz_high = star_properties->AGB_n_z - 1;
+ if (*iz_high >= star_properties->AGB_n_z)
+ *iz_high = star_properties->AGB_n_z - 1;
if (log_metallicity >= star_properties->yield_AGB.metallicity[0] &&
- log_metallicity <= star_properties->yield_AGB.metallicity[star_properties->AGB_n_z - 1])
+ log_metallicity <= star_properties->yield_AGB
+ .metallicity[star_properties->AGB_n_z - 1])
*dz = log_metallicity - star_properties->yield_AGB.metallicity[*iz_low];
else
*dz = 0;
- float deltaz = star_properties->yield_AGB.metallicity[*iz_high] - star_properties->yield_AGB.metallicity[*iz_low];
+ float deltaz = star_properties->yield_AGB.metallicity[*iz_high] -
+ star_properties->yield_AGB.metallicity[*iz_low];
if (deltaz > 0)
*dz = *dz / deltaz;
@@ -448,26 +489,32 @@ inline static void determine_bin_yield_AGB(int *iz_low, int *iz_high, float *dz,
}
}
-inline static void determine_bin_yield_SNII(int *iz_low, int *iz_high, float *dz, float log_metallicity,
- const struct stars_props *restrict star_properties){
+inline static void determine_bin_yield_SNII(
+ int* iz_low, int* iz_high, float* dz, float log_metallicity,
+ const struct stars_props* restrict star_properties) {
if (log_metallicity > log_min_metallicity) {
int j;
- for (j = 0; j < star_properties->SNII_n_z - 1 &&
- log_metallicity > star_properties->yield_SNII.metallicity[j + 1];
- j++);
+ for (j = 0;
+ j < star_properties->SNII_n_z - 1 &&
+ log_metallicity > star_properties->yield_SNII.metallicity[j + 1];
+ j++)
+ ;
*iz_low = j;
*iz_high = *iz_low + 1;
- if (*iz_high >= star_properties->SNII_n_z) *iz_high = star_properties->SNII_n_z - 1;
+ if (*iz_high >= star_properties->SNII_n_z)
+ *iz_high = star_properties->SNII_n_z - 1;
if (log_metallicity >= star_properties->yield_SNII.metallicity[0] &&
- log_metallicity <= star_properties->yield_SNII.metallicity[star_properties->SNII_n_z - 1])
+ log_metallicity <= star_properties->yield_SNII
+ .metallicity[star_properties->SNII_n_z - 1])
*dz = log_metallicity - star_properties->yield_SNII.metallicity[*iz_low];
else
*dz = 0;
- float deltaz = star_properties->yield_SNII.metallicity[*iz_high] - star_properties->yield_SNII.metallicity[*iz_low];
+ float deltaz = star_properties->yield_SNII.metallicity[*iz_high] -
+ star_properties->yield_SNII.metallicity[*iz_low];
if (deltaz > 0)
*dz = *dz / deltaz;
@@ -480,12 +527,11 @@ inline static void determine_bin_yield_SNII(int *iz_low, int *iz_high, float *dz
}
}
-
inline static void evolve_SNIa(float log10_min_mass, float log10_max_mass,
- const struct stars_props *restrict stars,
- struct spart *restrict sp,
- const struct unit_system* us,
- float star_age_Gyr, float dt_Gyr){
+ const struct stars_props* restrict stars,
+ struct spart* restrict sp,
+ const struct unit_system* us, float star_age_Gyr,
+ float dt_Gyr) {
/* Check if we're outside the mass range for SNIa */
if (log10_min_mass >= log10_SNIa_max_mass_msun) return;
@@ -493,7 +539,9 @@ inline static void evolve_SNIa(float log10_min_mass, float log10_max_mass,
/* check integration limits */
if (log10_max_mass > log10_SNIa_max_mass_msun) {
log10_max_mass = log10_SNIa_max_mass_msun;
- float lifetime_Gyr = lifetime_in_Gyr(exp(M_LN10 * log10_SNIa_max_mass_msun), sp->chemistry_data.metal_mass_fraction_total, stars);
+ float lifetime_Gyr =
+ lifetime_in_Gyr(exp(M_LN10 * 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;
}
@@ -503,34 +551,44 @@ inline static void evolve_SNIa(float log10_min_mass, float log10_max_mass,
sp->to_distribute.num_SNIa =
stars->SNIa_efficiency *
(exp(-star_age_Gyr / stars->SNIa_timescale) -
- exp(-(star_age_Gyr + dt_Gyr) / stars->SNIa_timescale))
- * sp->mass_init / stars->const_solar_mass;
+ exp(-(star_age_Gyr + dt_Gyr) / stars->SNIa_timescale)) *
+ sp->mass_init / stars->const_solar_mass;
/* compute total mass released by SNIa */
- sp->to_distribute.mass += sp->to_distribute.num_SNIa * stars->yield_SNIa_total_metals_SPH * stars->const_solar_mass;
+ sp->to_distribute.mass += sp->to_distribute.num_SNIa *
+ stars->yield_SNIa_total_metals_SPH *
+ stars->const_solar_mass;
/* compute mass fractions of each metal */
for (int i = 0; i < chemistry_element_count; i++) {
- sp->to_distribute.chemistry_data.metal_mass_fraction[i] += sp->to_distribute.num_SNIa * stars->yield_SNIa_SPH[i] * stars->const_solar_mass / sp->to_distribute.mass;
+ sp->to_distribute.chemistry_data.metal_mass_fraction[i] +=
+ sp->to_distribute.num_SNIa * stars->yield_SNIa_SPH[i] *
+ stars->const_solar_mass / sp->to_distribute.mass;
}
// For diagnostics according to Richard
- //if (stars->SNIa_mass_transfer == 1) {
+ // if (stars->SNIa_mass_transfer == 1) {
// for (i = 0; i < chemistry_element_count; i++) {
- // sp->metals_released[i] += num_of_SNIa_per_msun * stars->yield_SNIa_SPH[i];
+ // sp->metals_released[i] += num_of_SNIa_per_msun *
+ // stars->yield_SNIa_SPH[i];
// }
- // sp->chemistry_data.mass_from_SNIa += num_of_SNIa_per_msun * stars->yield_SNIa_total_metals_SPH;
- // sp->chemistry_data.metal_mass_fraction_from_SNIa += num_of_SNIa_per_msun * stars->yield_SNIa_total_metals_SPH;
+ // sp->chemistry_data.mass_from_SNIa += num_of_SNIa_per_msun *
+ // stars->yield_SNIa_total_metals_SPH;
+ // sp->chemistry_data.metal_mass_fraction_from_SNIa += num_of_SNIa_per_msun *
+ // stars->yield_SNIa_total_metals_SPH;
- // sp->metal_mass_released += num_of_SNIa_per_msun * stars->yield_SNIa_total_metals_SPH;
+ // sp->metal_mass_released += num_of_SNIa_per_msun *
+ // stars->yield_SNIa_total_metals_SPH;
- // // Make sure chemistry_element_Fe corresponds to the iron_index used in EAGLE!!!
- // sp->chemistry_data.iron_mass_fraction_from_SNIa += num_of_SNIa_per_msun * stars->yield_SNIa_SPH[chemistry_element_Fe];
+ // // Make sure chemistry_element_Fe corresponds to the iron_index used in
+ // EAGLE!!! sp->chemistry_data.iron_mass_fraction_from_SNIa +=
+ // num_of_SNIa_per_msun * stars->yield_SNIa_SPH[chemistry_element_Fe];
// /* metal_mass_released is the yield of ALL metals, not just the
// 11 tabulated in the code. SNIa remnants inject no H or He
- // so chemistry_data.mass_from_SNIa == chemistry_data.metal_mass_fraction_from_SNIa */
+ // so chemistry_data.mass_from_SNIa ==
+ // chemistry_data.metal_mass_fraction_from_SNIa */
//} else {
// sp->chemistry_data.iron_mass_fraction_from_SNIa = 0;
@@ -540,13 +598,13 @@ inline static void evolve_SNIa(float log10_min_mass, float log10_max_mass,
}
inline static void evolve_SNII(float log10_min_mass, float log10_max_mass,
- const struct stars_props *restrict stars,
- struct spart *restrict sp){
+ const struct stars_props* restrict stars,
+ struct spart* restrict sp) {
// come up with more descriptive index names
int ilow, ihigh, imass, i = 0;
/* determine integration range: make sure all these stars actually become SN
- * of type II */
+ * of type II */
if (log10_min_mass < log10_SNII_min_mass_msun)
log10_min_mass = log10_SNII_min_mass_msun;
@@ -558,44 +616,61 @@ inline static void evolve_SNII(float log10_min_mass, float log10_max_mass,
/* determine which mass bins will contribute */
determine_imf_bins(log10_min_mass, log10_max_mass, &ilow, &ihigh, stars);
- sp->to_distribute.num_SNII = integrate_imf(log10_min_mass, log10_max_mass, 0.0, 0, stars->stellar_yield,stars);
+ sp->to_distribute.num_SNII = integrate_imf(
+ log10_min_mass, log10_max_mass, 0.0, 0, stars->stellar_yield, stars);
/* determine yield of these bins (not equally spaced bins) */
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);
+ determine_bin_yield_SNII(&iz_low, &iz_high, &dz,
+ log10(sp->chemistry_data.metal_mass_fraction_total),
+ stars);
/* compute stellar_yield as function of mass */
float metals[chemistry_element_count], mass;
for (i = 0; i < chemistry_element_count; i++) {
for (imass = ilow; imass < ihigh + 1; imass++) {
- low_index_3d = row_major_index_3d(iz_low,i,imass,stars->SNII_n_z,chemistry_element_count,stars->SNII_n_mass);
- high_index_3d = row_major_index_3d(iz_high,i,imass,stars->SNII_n_z,chemistry_element_count,stars->SNII_n_mass);
- low_index_2d = row_major_index_2d(iz_low,imass,stars->SNII_n_z,stars->SNII_n_mass);
- high_index_2d = row_major_index_2d(iz_high,imass,stars->SNII_n_z,stars->SNII_n_mass);
- /* yield_SNII.SPH refers to elements produced, yield_SNII.ejecta_SPH refers
- * to elements already in star */
+ low_index_3d =
+ row_major_index_3d(iz_low, i, imass, stars->SNII_n_z,
+ chemistry_element_count, stars->SNII_n_mass);
+ high_index_3d =
+ row_major_index_3d(iz_high, i, imass, stars->SNII_n_z,
+ chemistry_element_count, stars->SNII_n_mass);
+ low_index_2d = row_major_index_2d(iz_low, imass, stars->SNII_n_z,
+ stars->SNII_n_mass);
+ high_index_2d = row_major_index_2d(iz_high, imass, stars->SNII_n_z,
+ stars->SNII_n_mass);
+ /* yield_SNII.SPH refers to elements produced, yield_SNII.ejecta_SPH
+ * refers to elements already in star */
stars->stellar_yield[imass] =
(1 - dz) * (stars->yield_SNII.SPH[low_index_3d] +
- sp->chemistry_data.metal_mass_fraction[i] * stars->yield_SNII.ejecta_SPH[low_index_2d]) +
+ sp->chemistry_data.metal_mass_fraction[i] *
+ stars->yield_SNII.ejecta_SPH[low_index_2d]) +
dz * (stars->yield_SNII.SPH[high_index_3d] +
- sp->chemistry_data.metal_mass_fraction[i] * stars->yield_SNII.ejecta_SPH[high_index_2d]);
+ sp->chemistry_data.metal_mass_fraction[i] *
+ stars->yield_SNII.ejecta_SPH[high_index_2d]);
}
- metals[i] = integrate_imf(log10_min_mass, log10_max_mass, 0.0, 2, stars->stellar_yield,stars);
+ metals[i] = integrate_imf(log10_min_mass, log10_max_mass, 0.0, 2,
+ stars->stellar_yield, stars);
}
for (imass = ilow; imass < ihigh + 1; imass++) {
- low_index_2d = row_major_index_2d(iz_low,imass,stars->SNII_n_z,stars->SNII_n_mass);
- high_index_2d = row_major_index_2d(iz_high,imass,stars->SNII_n_z,stars->SNII_n_mass);
+ low_index_2d =
+ row_major_index_2d(iz_low, imass, stars->SNII_n_z, stars->SNII_n_mass);
+ high_index_2d =
+ row_major_index_2d(iz_high, imass, stars->SNII_n_z, stars->SNII_n_mass);
stars->stellar_yield[imass] =
(1 - dz) * (stars->yield_SNII.total_metals_SPH[low_index_2d] +
- sp->chemistry_data.metal_mass_fraction_total * stars->yield_SNII.ejecta_SPH[low_index_2d]) +
+ sp->chemistry_data.metal_mass_fraction_total *
+ stars->yield_SNII.ejecta_SPH[low_index_2d]) +
dz * (stars->yield_SNII.total_metals_SPH[high_index_2d] +
- sp->chemistry_data.metal_mass_fraction_total * stars->yield_SNII.ejecta_SPH[high_index_2d]);
+ sp->chemistry_data.metal_mass_fraction_total *
+ stars->yield_SNII.ejecta_SPH[high_index_2d]);
}
- mass = integrate_imf(log10_min_mass, log10_max_mass, 0.0, 2, stars->stellar_yield, stars);
+ mass = integrate_imf(log10_min_mass, log10_max_mass, 0.0, 2,
+ stars->stellar_yield, stars);
/* yield normalization */
float norm0, norm1;
@@ -608,18 +683,23 @@ inline static void evolve_SNII(float log10_min_mass, float log10_max_mass,
/* get the total mass ejected from the table */
for (imass = ilow; imass < ihigh + 1; imass++) {
- low_index_2d = row_major_index_2d(iz_low,imass,stars->SNII_n_z,stars->SNII_n_mass);
- high_index_2d = row_major_index_2d(iz_high,imass,stars->SNII_n_z,stars->SNII_n_mass);
- stars->stellar_yield[imass] = (1 - dz) * stars->yield_SNII.ejecta_SPH[low_index_2d] +
- dz * stars->yield_SNII.ejecta_SPH[high_index_2d];
+ low_index_2d =
+ row_major_index_2d(iz_low, imass, stars->SNII_n_z, stars->SNII_n_mass);
+ high_index_2d =
+ row_major_index_2d(iz_high, imass, stars->SNII_n_z, stars->SNII_n_mass);
+ stars->stellar_yield[imass] =
+ (1 - dz) * stars->yield_SNII.ejecta_SPH[low_index_2d] +
+ dz * stars->yield_SNII.ejecta_SPH[high_index_2d];
}
- norm0 = integrate_imf(log10_min_mass, log10_max_mass, 0.0, 2, stars->stellar_yield, stars);
+ norm0 = integrate_imf(log10_min_mass, log10_max_mass, 0.0, 2,
+ stars->stellar_yield, stars);
/* compute the total mass ejected */
norm1 = mass + metals[chemistry_element_H] + metals[chemistry_element_He];
- /* Set normalisation factor. Note additional multiplication by the stellar initial mass as tables are per initial mass */
+ /* Set normalisation factor. Note additional multiplication by the stellar
+ * initial mass as tables are per initial mass */
const float norm_factor = norm0 / norm1 * sp->mass_init;
/* normalize the yields */
@@ -630,7 +710,11 @@ inline static void evolve_SNII(float log10_min_mass, float log10_max_mass,
sp->chemistry_data.mass_from_SNII += sp->metals_released[i];
}
sp->to_distribute.mass += mass * norm_factor;
- message("SNII mass to distribute %.5e initial mass %.5e norm mass %.5e mass %.5e norm_factor %.5e ", sp->to_distribute.mass, sp->mass_init, mass*norm_factor, mass, norm_factor);
+ message(
+ "SNII mass to distribute %.5e initial mass %.5e norm mass %.5e mass "
+ "%.5e norm_factor %.5e ",
+ sp->to_distribute.mass, sp->mass_init, mass * norm_factor, mass,
+ norm_factor);
sp->metal_mass_released += mass * norm_factor;
sp->chemistry_data.metal_mass_fraction_from_SNII += mass * norm_factor;
} else {
@@ -643,8 +727,8 @@ inline static void evolve_SNII(float log10_min_mass, float log10_max_mass,
}
inline static void evolve_AGB(float log10_min_mass, float log10_max_mass,
- const struct stars_props *restrict stars,
- struct spart *restrict sp){
+ const struct stars_props* restrict stars,
+ struct spart* restrict sp) {
// come up with more descriptive index names
int ilow, ihigh, imass, i = 0;
@@ -660,39 +744,55 @@ inline static void evolve_AGB(float log10_min_mass, float log10_max_mass,
/* determine yield of these bins (not equally spaced bins) */
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);
+ determine_bin_yield_AGB(&iz_low, &iz_high, &dz,
+ log10(sp->chemistry_data.metal_mass_fraction_total),
+ stars);
/* compute stellar_yield as function of mass */
float metals[chemistry_element_count], mass;
for (i = 0; i < chemistry_element_count; i++) {
for (imass = ilow; imass < ihigh + 1; imass++) {
- low_index_3d = row_major_index_3d(iz_low,i,imass,stars->AGB_n_z,chemistry_element_count,stars->AGB_n_mass);
- high_index_3d = row_major_index_3d(iz_high,i,imass,stars->AGB_n_z,chemistry_element_count,stars->AGB_n_mass);
- low_index_2d = row_major_index_2d(iz_low,imass,stars->AGB_n_z,stars->AGB_n_mass);
- high_index_2d = row_major_index_2d(iz_high,imass,stars->AGB_n_z,stars->AGB_n_mass);
+ low_index_3d =
+ row_major_index_3d(iz_low, i, imass, stars->AGB_n_z,
+ chemistry_element_count, stars->AGB_n_mass);
+ high_index_3d =
+ row_major_index_3d(iz_high, i, imass, stars->AGB_n_z,
+ chemistry_element_count, stars->AGB_n_mass);
+ low_index_2d =
+ row_major_index_2d(iz_low, imass, stars->AGB_n_z, stars->AGB_n_mass);
+ high_index_2d =
+ row_major_index_2d(iz_high, imass, stars->AGB_n_z, stars->AGB_n_mass);
/* yield_AGB.SPH refers to elements produced, yield_AGB.ejecta_SPH refers
* to elements already in star */
stars->stellar_yield[imass] =
(1 - dz) * (stars->yield_AGB.SPH[low_index_3d] +
- sp->chemistry_data.metal_mass_fraction[i] * stars->yield_AGB.ejecta_SPH[low_index_2d]) +
+ sp->chemistry_data.metal_mass_fraction[i] *
+ stars->yield_AGB.ejecta_SPH[low_index_2d]) +
dz * (stars->yield_AGB.SPH[high_index_3d] +
- sp->chemistry_data.metal_mass_fraction[i] * stars->yield_AGB.ejecta_SPH[high_index_2d]);
+ sp->chemistry_data.metal_mass_fraction[i] *
+ stars->yield_AGB.ejecta_SPH[high_index_2d]);
}
- metals[i] = integrate_imf(log10_min_mass, log10_max_mass, 0.0, 2, stars->stellar_yield,stars);
+ metals[i] = integrate_imf(log10_min_mass, log10_max_mass, 0.0, 2,
+ stars->stellar_yield, stars);
}
for (imass = ilow; imass < ihigh + 1; imass++) {
- low_index_2d = row_major_index_2d(iz_low,imass,stars->AGB_n_z,stars->AGB_n_mass);
- high_index_2d = row_major_index_2d(iz_high,imass,stars->AGB_n_z,stars->AGB_n_mass);
+ low_index_2d =
+ row_major_index_2d(iz_low, imass, stars->AGB_n_z, stars->AGB_n_mass);
+ high_index_2d =
+ row_major_index_2d(iz_high, imass, stars->AGB_n_z, stars->AGB_n_mass);
stars->stellar_yield[imass] =
(1 - dz) * (stars->yield_AGB.total_metals_SPH[low_index_2d] +
- sp->chemistry_data.metal_mass_fraction_total * stars->yield_AGB.ejecta_SPH[low_index_2d]) +
+ sp->chemistry_data.metal_mass_fraction_total *
+ stars->yield_AGB.ejecta_SPH[low_index_2d]) +
dz * (stars->yield_AGB.total_metals_SPH[high_index_2d] +
- sp->chemistry_data.metal_mass_fraction_total * stars->yield_AGB.ejecta_SPH[high_index_2d]);
+ sp->chemistry_data.metal_mass_fraction_total *
+ stars->yield_AGB.ejecta_SPH[high_index_2d]);
}
- mass = integrate_imf(log10_min_mass, log10_max_mass, 0.0, 2, stars->stellar_yield, stars);
+ mass = integrate_imf(log10_min_mass, log10_max_mass, 0.0, 2,
+ stars->stellar_yield, stars);
/* yield normalization */
float norm0, norm1;
@@ -705,30 +805,42 @@ inline static void evolve_AGB(float log10_min_mass, float log10_max_mass,
/* get the total mass ejected from the table */
for (imass = ilow; imass < ihigh + 1; imass++) {
- low_index_2d = row_major_index_2d(iz_low,imass,stars->AGB_n_z,stars->AGB_n_mass);
- high_index_2d = row_major_index_2d(iz_high,imass,stars->AGB_n_z,stars->AGB_n_mass);
- stars->stellar_yield[imass] = (1 - dz) * stars->yield_AGB.ejecta_SPH[low_index_2d] +
- dz * stars->yield_AGB.ejecta_SPH[high_index_2d];
+ low_index_2d =
+ row_major_index_2d(iz_low, imass, stars->AGB_n_z, stars->AGB_n_mass);
+ high_index_2d =
+ row_major_index_2d(iz_high, imass, stars->AGB_n_z, stars->AGB_n_mass);
+ stars->stellar_yield[imass] =
+ (1 - dz) * stars->yield_AGB.ejecta_SPH[low_index_2d] +
+ dz * stars->yield_AGB.ejecta_SPH[high_index_2d];
}
- norm0 = integrate_imf(log10_min_mass, log10_max_mass, 0.0, 2, stars->stellar_yield, stars);
+ norm0 = integrate_imf(log10_min_mass, log10_max_mass, 0.0, 2,
+ stars->stellar_yield, stars);
/* compute the total mass ejected */
norm1 = mass + metals[chemistry_element_H] + metals[chemistry_element_He];
- /* Set normalisation factor. Note additional multiplication by the stellar initial mass as tables are per initial mass */
- const float norm_factor = norm0/norm1 * sp->mass_init;
+ /* Set normalisation factor. Note additional multiplication by the stellar
+ * initial mass as tables are per initial mass */
+ const float norm_factor = norm0 / norm1 * sp->mass_init;
/* normalize the yields (Copied from SNII) */
if (norm1 > 0) {
for (i = 0; i < chemistry_element_count; i++) {
sp->metals_released[i] += metals[i] * norm_factor;
- // This increment is coppied from EAGLE, however note that it is different from SNII case. Investigate?
+ // This increment is coppied from EAGLE, however note that it is different
+ // from SNII case. Investigate?
sp->chemistry_data.mass_from_AGB += metals[i] * norm_factor;
}
- sp->to_distribute.mass += (mass + metals[chemistry_element_H] + metals[chemistry_element_He]) * norm_factor;
- message("AGB mass to distribute %.5e initial mass %.5e norm mass %.5e mass %.5e norm_factor %.5e ", sp->to_distribute.mass, sp->mass_init, mass*norm_factor, mass, norm_factor);
+ sp->to_distribute.mass +=
+ (mass + metals[chemistry_element_H] + metals[chemistry_element_He]) *
+ norm_factor;
+ message(
+ "AGB mass to distribute %.5e initial mass %.5e norm mass %.5e mass "
+ "%.5e norm_factor %.5e ",
+ sp->to_distribute.mass, sp->mass_init, mass * norm_factor, mass,
+ norm_factor);
sp->metal_mass_released += mass * norm_factor;
sp->chemistry_data.metal_mass_fraction_from_AGB += mass * norm_factor;
} else {
@@ -737,43 +849,55 @@ inline static void evolve_AGB(float log10_min_mass, float log10_max_mass,
}
// Analogue of eagle_do_stellar_evolution...
-inline static void compute_stellar_evolution(const struct stars_props *restrict star_properties,
- struct spart *restrict sp,
- const struct unit_system* us,
- float age, float dt) {
-
+inline static void compute_stellar_evolution(
+ const struct stars_props* restrict star_properties,
+ struct spart* restrict sp, const struct unit_system* us, float age,
+ float dt) {
+
// Convert dt from internal units to Gyr.
const float Gyr_in_cgs = 3.154e16;
- float dt_Gyr = dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME) / Gyr_in_cgs;
- float star_age_Gyr = age * units_cgs_conversion_factor(us, UNIT_CONV_TIME) / Gyr_in_cgs; // This should be same as age_of_star_in_Gyr_begstep in EAGLE, check. Also, make sure it works with cosmology, might need to use birth_scale_factor.
+ float dt_Gyr =
+ dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME) / Gyr_in_cgs;
+ float star_age_Gyr =
+ age * units_cgs_conversion_factor(us, UNIT_CONV_TIME) /
+ Gyr_in_cgs; // This should be same as age_of_star_in_Gyr_begstep in
+ // EAGLE, check. Also, make sure it works with cosmology,
+ // might need to use birth_scale_factor.
// set max and min mass of dying stars
- float log10_max_dying_mass_msun =
- log10(dying_mass_msun(star_age_Gyr, sp->chemistry_data.metal_mass_fraction_total, star_properties));
- float log10_min_dying_mass_msun = log10(
- dying_mass_msun(star_age_Gyr + dt_Gyr, sp->chemistry_data.metal_mass_fraction_total, star_properties));
+ float log10_max_dying_mass_msun = log10(dying_mass_msun(
+ star_age_Gyr, sp->chemistry_data.metal_mass_fraction_total,
+ star_properties));
+ float log10_min_dying_mass_msun = log10(dying_mass_msun(
+ star_age_Gyr + dt_Gyr, sp->chemistry_data.metal_mass_fraction_total,
+ star_properties));
- if (log10_min_dying_mass_msun > log10_max_dying_mass_msun) error("min dying mass is greater than max dying mass");
+ 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 */
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,us,star_age_Gyr,dt_Gyr);
- //evolve_SNII(log10_min_dying_mass_msun,log10_max_dying_mass_msun,star_properties,sp);
- evolve_AGB(log10_min_dying_mass_msun,log10_max_dying_mass_msun,star_properties,sp);
-
- sp->to_distribute.chemistry_data.metal_mass_fraction_total = 1.f - sp->to_distribute.chemistry_data.metal_mass_fraction[0] - sp->to_distribute.chemistry_data.metal_mass_fraction[1];
-
+ // evolve_SNIa(log10_min_dying_mass_msun,log10_max_dying_mass_msun,star_properties,sp,us,star_age_Gyr,dt_Gyr);
+ // evolve_SNII(log10_min_dying_mass_msun,log10_max_dying_mass_msun,star_properties,sp);
+ evolve_AGB(log10_min_dying_mass_msun, log10_max_dying_mass_msun,
+ star_properties, sp);
+
+ sp->to_distribute.chemistry_data.metal_mass_fraction_total =
+ 1.f - sp->to_distribute.chemistry_data.metal_mass_fraction[0] -
+ sp->to_distribute.chemistry_data.metal_mass_fraction[1];
}
inline static float compute_SNe(struct spart* sp,
- const struct stars_props* stars_properties,
- float age, double dt) {
- if (age <= stars_properties->SNII_wind_delay && age + dt > stars_properties->SNII_wind_delay) {
+ const struct stars_props* stars_properties,
+ float age, double dt) {
+ if (age <= stars_properties->SNII_wind_delay &&
+ age + dt > stars_properties->SNII_wind_delay) {
// ALEXEI: commented for debugging
- //return stars_properties->num_SNII_per_msun * sp->mass_init / stars_properties->const_solar_mass;
+ // return stars_properties->num_SNII_per_msun * sp->mass_init /
+ // stars_properties->const_solar_mass;
return 0;
} else {
return 0;
@@ -792,37 +916,38 @@ inline static float compute_SNe(struct spart* sp,
*/
__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) {
-
- // Set birth time for testing purposes
- sp->birth_time = 0;
- float star_age = current_time - sp->birth_time;
-
- sp->to_distribute.num_SNIa = 0;
- sp->to_distribute.num_SNII = 0;
- sp->to_distribute.mass = 0;
-
- // Set elements released to zero
- for(int i = 0; i < chemistry_element_count; i++) sp->metals_released[i] = 0;
- sp->metal_mass_released = 0;
- sp->chemistry_data.mass_from_AGB = 0;
- sp->chemistry_data.metal_mass_fraction_from_AGB = 0;
- sp->chemistry_data.mass_from_SNII = 0;
- sp->chemistry_data.metal_mass_fraction_from_SNII = 0;
- sp->chemistry_data.mass_from_SNIa = 0;
- sp->chemistry_data.metal_mass_fraction_from_SNIa = 0;
- sp->chemistry_data.iron_mass_fraction_from_SNIa = 0;
-
- // Evolve the star
- compute_stellar_evolution(stars_properties, sp, us, star_age, dt);
-
- /* Compute the number of type II SNe that went off */
- sp->to_distribute.num_SNe = compute_SNe(sp,stars_properties,star_age,dt);
-
+ const struct cosmology* cosmo, const struct unit_system* us,
+ float current_time, double dt) {
+
+ // Set birth time for testing purposes
+ sp->birth_time = 0;
+ float star_age = current_time - sp->birth_time;
+
+ sp->to_distribute.num_SNIa = 0;
+ sp->to_distribute.num_SNII = 0;
+ sp->to_distribute.mass = 0;
+
+ // Set elements released to zero
+ for (int i = 0; i < chemistry_element_count; i++) sp->metals_released[i] = 0;
+ sp->metal_mass_released = 0;
+ sp->chemistry_data.mass_from_AGB = 0;
+ sp->chemistry_data.metal_mass_fraction_from_AGB = 0;
+ sp->chemistry_data.mass_from_SNII = 0;
+ sp->chemistry_data.metal_mass_fraction_from_SNII = 0;
+ sp->chemistry_data.mass_from_SNIa = 0;
+ sp->chemistry_data.metal_mass_fraction_from_SNIa = 0;
+ sp->chemistry_data.iron_mass_fraction_from_SNIa = 0;
+
+ // Evolve the star
+ compute_stellar_evolution(stars_properties, sp, us, star_age, dt);
+
+ /* Compute the number of type II SNe that went off */
+ sp->to_distribute.num_SNe = compute_SNe(sp, stars_properties, star_age, dt);
}
-inline static void stars_evolve_init(struct swift_params *params, struct stars_props* restrict stars){
-
+inline static void stars_evolve_init(struct swift_params* params,
+ struct stars_props* restrict stars) {
+
stars->SNIa_n_elements = 42;
stars->SNII_n_mass = 11;
stars->SNII_n_elements = 11;
@@ -834,19 +959,21 @@ inline static void stars_evolve_init(struct swift_params *params, struct stars_p
stars->lifetimes.n_z = 6;
stars->element_name_length = 15;
- /* Turn on AGB and SNII mass transfer (Do we really need this?
- * Should these maybe always be on? If not on they effectively
+ /* Turn on AGB and SNII mass transfer (Do we really need this?
+ * Should these maybe always be on? If not on they effectively
* turn off SNII and AGB evolution.) */
stars->AGB_mass_transfer = 1;
stars->SNII_mass_transfer = 1;
/* Yield table filepath */
- parser_get_param_string(params, "EagleStellarEvolution:filename", stars->yield_table_path);
- parser_get_param_string(params, "EagleStellarEvolution:imf_model", stars->IMF_Model);
+ parser_get_param_string(params, "EagleStellarEvolution:filename",
+ stars->yield_table_path);
+ parser_get_param_string(params, "EagleStellarEvolution:imf_model",
+ stars->IMF_Model);
/* Allocate yield tables */
allocate_yield_tables(stars);
-
+
/* Set factors for each element adjusting SNII yield */
stars->typeII_factor[0] = 1.f;
stars->typeII_factor[1] = 1.f;
@@ -867,21 +994,24 @@ inline static void stars_evolve_init(struct swift_params *params, struct stars_p
/* Set yield_mass_bins array */
const float lm_min = log10(imf_min_mass_msun); /* min mass in solar masses */
const float lm_max = log10(imf_max_mass_msun); /* max mass in solar masses */
- const float dlm = (lm_max - lm_min) / (n_mass_bins - 1);
+ const float dlm = (lm_max - lm_min) / (n_mass_bins - 1);
// ALEXEI: does yield_mass_bins really have to be double?
- for (int i = 0; i < n_mass_bins; i++) stars->yield_mass_bins[i] = dlm * i + lm_min;
+ for (int i = 0; i < n_mass_bins; i++)
+ stars->yield_mass_bins[i] = dlm * i + lm_min;
- /* Further calculation on tables to convert them to log10 and compute yields for each element */
+ /* Further calculation on tables to convert them to log10 and compute yields
+ * for each element */
compute_yields(stars);
-
+
/* Further calculation on tables to compute ejecta */
compute_ejecta(stars);
/* Calculate number of type II SN per solar mass */
- stars->num_SNII_per_msun = integrate_imf(stars->log10_SNII_min_mass_msun,stars->log10_SNII_max_mass_msun, 0, 0, stars->stellar_yield, stars);
+ stars->num_SNII_per_msun = integrate_imf(stars->log10_SNII_min_mass_msun,
+ stars->log10_SNII_max_mass_msun, 0,
+ 0, stars->stellar_yield, stars);
message("initialized stellar feedback");
-
}
/**
diff --git a/src/stars/EAGLE/stars_iact.h b/src/stars/EAGLE/stars_iact.h
index ea8c3edac9feb8216d78f8d01a59e08068ce4e82..c9f13d1f02ca055e11c08e82c33ca43a9656f9c2 100644
--- a/src/stars/EAGLE/stars_iact.h
+++ b/src/stars/EAGLE/stars_iact.h
@@ -36,12 +36,11 @@
* @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) {
+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) {
/* Get the gas mass. */
const float mj = hydro_get_mass(pj);
@@ -56,7 +55,7 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
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;
@@ -67,11 +66,12 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
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);// / stars_properties->const_solar_mass;
+ si->ngb_mass += hydro_get_mass(pj); // / stars_properties->const_solar_mass;
+
+ /* Add contribution of pj to normalisation of kernel (TODO: IMPROVE COMMENT?)
+ */
+ si->omega_normalisation_inv += wj / hydro_get_physical_density(pj, cosmo);
- /* Add contribution of pj to normalisation of kernel (TODO: IMPROVE COMMENT?) */
- si->omega_normalisation_inv += wj / hydro_get_physical_density(pj,cosmo);
-
/* Compute contribution to the density */
si->rho_gas += mj * wi;
@@ -84,7 +84,7 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
}
/**
- * @brief Increases thermal energy of particle due
+ * @brief Increases thermal energy of particle due
* to feedback by specified amount
*
* @param du change in internal energy
@@ -92,15 +92,16 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
* @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) {
+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.
+ // 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;
}
@@ -117,15 +118,15 @@ static inline void thermal_feedback(float du, struct part * restrict p,
* @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
+ * @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) {
+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. */
@@ -138,63 +139,90 @@ runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
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);
+ 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;
- hydro_set_mass(pj,new_mass);
+ float new_mass = current_mass + si->to_distribute.mass * omega_frac;
+ hydro_set_mass(pj, new_mass);
/* Decrease the mass of star particle */
- si->mass -= si->to_distribute.mass*omega_frac;
+ 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;
-
+ 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;
+ 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;
+ 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;
+ 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;
+ 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;
+ 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;
+ 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;
+ 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;
+ 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++) {
@@ -204,32 +232,44 @@ runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
/* Energy feedback */
float heating_probability = -1.f, du = 0.f, d_energy = 0.f;
- d_energy = si->to_distribute.mass * (stars_properties->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);
+ d_energy =
+ si->to_distribute.mass *
+ (stars_properties->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 (stars_properties->continuous_heating) {
- // We're doing ONLY continuous heating
- d_energy += si->to_distribute.num_SNIa * stars_properties->total_energy_SNe * omega_frac * si->mass_init;
- du = d_energy/hydro_get_mass(pj);
- thermal_feedback(du,pj,xp,cosmo);
+ // We're doing ONLY continuous heating
+ d_energy += si->to_distribute.num_SNIa *
+ stars_properties->total_energy_SNe * omega_frac * si->mass_init;
+ du = d_energy / hydro_get_mass(pj);
+ thermal_feedback(du, pj, xp, cosmo);
} else {
// We're doing stochastic heating
- heating_probability = stars_properties->total_energy_SNe / stars_properties->temp_to_u_factor * si->to_distribute.num_SNe *
+ heating_probability = stars_properties->total_energy_SNe /
+ stars_properties->temp_to_u_factor *
+ si->to_distribute.num_SNe *
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;
+ 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;
+ 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);
+ 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 du/ini %.5e", pj->id, heating_probability, random_num, du, du/hydro_get_physical_internal_energy(pj,xp,cosmo));
+ message(
+ "we did some heating! id %llu probability %.5e random_num %.5e du %.5e "
+ "du/ini %.5e",
+ pj->id, heating_probability, random_num, du,
+ du / hydro_get_physical_internal_energy(pj, xp, cosmo));
thermal_feedback(du, pj, xp, cosmo);
}
-
}
#endif /* SWIFT_EAGLE_STARS_IACT_H */
diff --git a/src/stars/EAGLE/stars_io.h b/src/stars/EAGLE/stars_io.h
index 3819d4dafb48a676d5f9c55ec5cd5a62b3620452..3614370af80e4852a2433c10171c51eb302abfc7 100644
--- a/src/stars/EAGLE/stars_io.h
+++ b/src/stars/EAGLE/stars_io.h
@@ -103,7 +103,7 @@ INLINE static void stars_props_init(struct stars_props *sp,
const struct unit_system *us,
struct swift_params *params,
const struct hydro_props *p,
- const struct cosmology *cosmo) {
+ const struct cosmology *cosmo) {
/* Kernel properties */
sp->eta_neighbours = parser_get_opt_param_float(
@@ -135,38 +135,46 @@ INLINE static void stars_props_init(struct stars_props *sp,
else
sp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
- sp->stellar_lifetime_flag = parser_get_opt_param_int(
- params, "EAGLEFeedback:lifetime_flag", 0);
+ sp->stellar_lifetime_flag =
+ parser_get_opt_param_int(params, "EAGLEFeedback:lifetime_flag", 0);
- sp->SNIa_timescale = parser_get_opt_param_float(
- params, "EAGLEFeedback:SNIa_timescale", 2.f);
- sp->SNIa_efficiency = parser_get_opt_param_float(
+ sp->SNIa_timescale =
+ parser_get_opt_param_float(params, "EAGLEFeedback:SNIa_timescale", 2.f);
+ sp->SNIa_efficiency = parser_get_opt_param_float(
params, "EAGLEFeedback:SNIa_efficiency", 2.e-3);
- sp->continuous_heating = parser_get_opt_param_int(
+ sp->continuous_heating = parser_get_opt_param_int(
params, "EAGLEFeedback:continuous_heating_switch", 0);
const float Gyr_in_cgs = 3.154e16;
- sp->SNII_wind_delay = parser_get_opt_param_float(
- params, "EAGLEFeedback:SNII_wind_delay_Gyr", 0.03) * Gyr_in_cgs / units_cgs_conversion_factor(us,UNIT_CONV_TIME);
+ sp->SNII_wind_delay = parser_get_opt_param_float(
+ params, "EAGLEFeedback:SNII_wind_delay_Gyr", 0.03) *
+ Gyr_in_cgs /
+ units_cgs_conversion_factor(us, UNIT_CONV_TIME);
// ALEXEI: find out where this gets set in EAGLE
sp->SNIa_energy_fraction = 1.0;
/* Set the temperature to use in stochastic heating */
- sp->SNe_deltaT_desired = 3.16228e7 / units_cgs_conversion_factor(us,UNIT_CONV_TEMPERATURE);
+ sp->SNe_deltaT_desired =
+ 3.16228e7 / 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
+ const float ejecta_velocity =
+ 1.0e6 / units_cgs_conversion_factor(
+ us, UNIT_CONV_SPEED); // EAGLE parameter is 10 km/s
sp->ejecta_specific_thermal_energy = 0.5 * ejecta_velocity * ejecta_velocity;
/* Energy released by supernova */
- sp->total_energy_SNe = 1.0e51/units_cgs_conversion_factor(us,UNIT_CONV_ENERGY);
+ sp->total_energy_SNe =
+ 1.0e51 / units_cgs_conversion_factor(us, UNIT_CONV_ENERGY);
/* 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);
+ sp->temp_to_u_factor =
+ phys_const->const_boltzmann_k /
+ (p->mu_ionised * (hydro_gamma_minus_one)*phys_const->const_proton_mass);
/* Calculate number of type II SN per solar mass */
- sp->log10_SNII_min_mass_msun = 0.77815125f; // log10(6).
- sp->log10_SNII_max_mass_msun = 2.f; // log10(100).
+ sp->log10_SNII_min_mass_msun = 0.77815125f; // log10(6).
+ sp->log10_SNII_max_mass_msun = 2.f; // log10(100).
/* Copy over solar mass */
sp->const_solar_mass = phys_const->const_solar_mass;
diff --git a/src/stars/EAGLE/stars_part.h b/src/stars/EAGLE/stars_part.h
index e82fbb94ab3fed65720faa0ea85789c55ba2ca64..3acc2bf6f7c9d7dc54f5c9868aa526bd22176d29 100644
--- a/src/stars/EAGLE/stars_part.h
+++ b/src/stars/EAGLE/stars_part.h
@@ -38,14 +38,14 @@ struct spart {
long long id;
/*! Pointer to corresponding gravity part. */
- struct gpart* gpart;
+ 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];
@@ -102,7 +102,8 @@ struct spart {
} to_distribute;
- /* kernel normalisation factor (equivalent to metalweight_norm in eagle_enrich.c:811, TODO: IMPROVE COMMENT) */
+ /* 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 */
@@ -176,19 +177,18 @@ struct yield_table {
struct lifetime_table {
/* number of elements, mass, and initial metallicity bins */
int n_mass;
- int n_z;
+ int n_z;
/* table of masses */
- double *mass;
+ double *mass;
/* table of metallicities */
- double *metallicity;
+ double *metallicity;
/* table of lifetimes depending on mass an metallicity */
- double **dyingtime;
+ double **dyingtime;
};
-
/**
* @brief Contains all the constants and parameters of the stars scheme
*/
@@ -215,7 +215,8 @@ struct stars_props {
/* 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) */
+ /* 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 */
@@ -233,7 +234,8 @@ struct stars_props {
/* 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) */
+ /* Number of supernovae per solar mass (used only for testing in const
+ * feedback model) */
float sn_per_msun;
/* Solar mass */
@@ -256,12 +258,12 @@ struct stars_props {
double *yield_SNIa_SPH;
double yield_SNIa_total_metals_SPH;
double *yields_SNIa;
-
+
/* Parameters to SNIa enrichment model */
int SNIa_mode;
float SNIa_efficiency;
float SNIa_timescale;
-
+
/* Mass transfer due to enrichment */
int SNIa_mass_transfer;
int SNII_mass_transfer;
@@ -299,12 +301,14 @@ struct stars_props {
struct lifetime_table lifetimes;
/* Flag defining stellar lifetime model */
- int stellar_lifetime_flag; // 0 for padovani & matteucci 1993, 1 for maeder & meynet 1989, 2 for portinari et al. 1998.
+ int stellar_lifetime_flag; // 0 for padovani & matteucci 1993, 1 for maeder &
+ // meynet 1989, 2 for portinari et al. 1998.
/* Location of yield tables */
char yield_table_path[50];
- /* Array for storing stellar yields for computation in IMF (clarify name, comment) */
+ /* Array for storing stellar yields for computation in IMF (clarify name,
+ * comment) */
float *stellar_yield;
/* number of type II supernovae per solar mass */
@@ -316,7 +320,6 @@ struct stars_props {
/* log10 min, max mass in SNII tables */
float log10_SNII_min_mass_msun;
float log10_SNII_max_mass_msun;
-
};
#endif /* SWIFT_EAGLE_STAR_PART_H */
diff --git a/src/stars/EAGLE/yield_tables.h b/src/stars/EAGLE/yield_tables.h
index 96910cce7ea509b3ff57949059bcbcd802c022af..6529cff297135497c262c83d8de947cf26fce4d6 100644
--- a/src/stars/EAGLE/yield_tables.h
+++ b/src/stars/EAGLE/yield_tables.h
@@ -19,13 +19,14 @@
#ifndef SWIFT_EAGLE_STARS_YIELD_TABLES_H
#define SWIFT_EAGLE_STARS_YIELD_TABLES_H
-#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
+#include <gsl/gsl_math.h>
#include <gsl/gsl_spline.h>
#include "chemistry.h"
static const float log_min_metallicity = -20;
-static const int n_mass_bins = 200; // temporary, put in correct value and move elsewhere.
+static const int n_mass_bins =
+ 200; // temporary, put in correct value and move elsewhere.
// Temporary, these functions need to be somewhere else
inline static char *mystrdup(const char *s) {
@@ -61,9 +62,11 @@ __attribute__((always_inline)) INLINE int row_major_index_3d(int i, int j,
return i * ny * nz + j * nz + k;
}
-inline static int get_element_index(const char *element_name, char **element_array, int n_elements){
-
- /* Compare element name we are trying to index to every name in element array */
+inline static int get_element_index(const char *element_name,
+ char **element_array, int n_elements) {
+
+ /* Compare element name we are trying to index to every name in element array
+ */
for (int i = 0; i < n_elements; i++) {
if (strcmp(element_array[i], element_name) == 0) return i;
}
@@ -72,13 +75,13 @@ inline static int get_element_index(const char *element_name, char **element_arr
return -1;
}
-inline static void read_yield_tables(struct stars_props *restrict stars){
+inline static void read_yield_tables(struct stars_props *restrict stars) {
#ifdef HAVE_HDF5
int i, j, k, index;
char fname[256], setname[100];
-
+
hid_t file_id, dataset, datatype;
herr_t status;
@@ -94,7 +97,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->SNIa_element_names);
+ stars->SNIa_element_names);
if (status < 0) error("error reading SNIa element names");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -102,19 +105,19 @@ inline static void read_yield_tables(struct stars_props *restrict stars){
if (status < 0) error("error closing datatype");
// What is this for? Copied from EAGLE
- //for (i = 0; i < stars->SNIa_n_elements; i++) {
+ // for (i = 0; i < stars->SNIa_n_elements; i++) {
// stars->SNIa_element_names[i] = mystrdup(stars->SNIa_element_names[i]);
//}
dataset = H5Dopen(file_id, "Yield", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->yields_SNIa);
+ stars->yields_SNIa);
if (status < 0) error("error reading SNIa yields");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
dataset = H5Dopen(file_id, "Total_Metals", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- &stars->yield_SNIa_total_metals_SPH);
+ &stars->yield_SNIa_total_metals_SPH);
if (status < 0) error("error reading SNIa total metal");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -135,7 +138,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->SNII_element_names);
+ stars->SNII_element_names);
if (status < 0) error("error reading SNII element names");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -143,19 +146,19 @@ inline static void read_yield_tables(struct stars_props *restrict stars){
if (status < 0) error("error closing datatype");
// What is this for (again)? copied from EAGLE
- //for (i = 0; i < stars->SNII_n_elements; i++)
+ // for (i = 0; i < stars->SNII_n_elements; i++)
// yieldsSNII.ElementName[i] = mystrdup(yieldsSNII.ElementName[i]);
dataset = H5Dopen(file_id, "Masses", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->yield_SNII.mass);
+ stars->yield_SNII.mass);
if (status < 0) error("error reading SNII masses");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
dataset = H5Dopen(file_id, "Metallicities", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->yield_SNII.metallicity);
+ stars->yield_SNII.metallicity);
if (status < 0) error("error reading SNII metallicities");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -180,14 +183,15 @@ inline static void read_yield_tables(struct stars_props *restrict stars){
sprintf(setname, "/Yields/%s/Yield", tempname1[i]);
dataset = H5Dopen(file_id, setname, H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- tempyield1);
+ tempyield1);
if (status < 0) error("error reading SNII yield");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
sprintf(setname, "/Yields/%s/Ejected_mass", tempname1[i]);
dataset = H5Dopen(file_id, setname, H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, tempej1);
+ status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+ tempej1);
if (status < 0) error("error reading SNII ejected masses");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -195,18 +199,19 @@ inline static void read_yield_tables(struct stars_props *restrict stars){
sprintf(setname, "/Yields/%s/Total_Metals", tempname1[i]);
dataset = H5Dopen(file_id, setname, H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- tempmet1);
+ tempmet1);
if (status < 0) error("error reading SNII total metals");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
for (k = 0; k < stars->SNII_n_mass; k++) {
- index = row_major_index_2d(i,k,stars->SNII_n_z,stars->SNII_n_mass);
+ index = row_major_index_2d(i, k, stars->SNII_n_z, stars->SNII_n_mass);
stars->yield_SNII.ejecta[index] = tempej1[k];
stars->yield_SNII.total_metals[index] = tempmet1[k];
for (j = 0; j < stars->SNII_n_elements; j++) {
- index = row_major_index_3d(i,j,k,stars->SNII_n_z,stars->SNII_n_elements,stars->SNII_n_mass);
+ index = row_major_index_3d(i, j, k, stars->SNII_n_z,
+ stars->SNII_n_elements, stars->SNII_n_mass);
stars->yield_SNII.yield[index] = tempyield1[j][k];
}
}
@@ -228,7 +233,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->AGB_element_names);
+ stars->AGB_element_names);
if (status < 0) error("error reading AGB element names");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -236,19 +241,19 @@ inline static void read_yield_tables(struct stars_props *restrict stars){
if (status < 0) error("error closing datatype");
// What is this for (again)? copied from EAGLE
- //for (i = 0; i < stars->AGB_n_elements; i++)
+ // for (i = 0; i < stars->AGB_n_elements; i++)
// yieldsAGB.ElementName[i] = mystrdup(yieldsAGB.ElementName[i]);
dataset = H5Dopen(file_id, "Masses", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->yield_AGB.mass);
+ stars->yield_AGB.mass);
if (status < 0) error("error reading AGB masses");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
dataset = H5Dopen(file_id, "Metallicities", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->yield_AGB.metallicity);
+ stars->yield_AGB.metallicity);
if (status < 0) error("error reading AGB metallicities");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -273,14 +278,15 @@ inline static void read_yield_tables(struct stars_props *restrict stars){
sprintf(setname, "/Yields/%s/Yield", tempname2[i]);
dataset = H5Dopen(file_id, setname, H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- tempyield2);
+ tempyield2);
if (status < 0) error("error reading AGB yield");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
sprintf(setname, "/Yields/%s/Ejected_mass", tempname2[i]);
dataset = H5Dopen(file_id, setname, H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, tempej2);
+ status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+ tempej2);
if (status < 0) error("error reading AGB ejected masses");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -288,18 +294,19 @@ inline static void read_yield_tables(struct stars_props *restrict stars){
sprintf(setname, "/Yields/%s/Total_Metals", tempname2[i]);
dataset = H5Dopen(file_id, setname, H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- tempmet2);
+ tempmet2);
if (status < 0) error("error reading AGB total metals");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
for (k = 0; k < stars->AGB_n_mass; k++) {
- index = row_major_index_2d(i,k,stars->AGB_n_z,stars->AGB_n_mass);
+ index = row_major_index_2d(i, k, stars->AGB_n_z, stars->AGB_n_mass);
stars->yield_AGB.ejecta[index] = tempej2[k];
stars->yield_AGB.total_metals[index] = tempmet2[k];
for (j = 0; j < stars->AGB_n_elements; j++) {
- index = row_major_index_3d(i,j,k,stars->AGB_n_z,stars->AGB_n_elements,stars->AGB_n_mass);
+ index = row_major_index_3d(i, j, k, stars->AGB_n_z,
+ stars->AGB_n_elements, stars->AGB_n_mass);
stars->yield_AGB.yield[index] = tempyield2[j][k];
}
}
@@ -318,14 +325,14 @@ inline static void read_yield_tables(struct stars_props *restrict stars){
dataset = H5Dopen(file_id, "Masses", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->lifetimes.mass);
+ stars->lifetimes.mass);
if (status < 0) error("error reading lifetime table masses");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
dataset = H5Dopen(file_id, "Metallicities", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- stars->lifetimes.metallicity);
+ stars->lifetimes.metallicity);
if (status < 0) error("error reading AGB metallicities");
status = H5Dclose(dataset);
if (status < 0) error("error closing dataset");
@@ -338,136 +345,186 @@ inline static void read_yield_tables(struct stars_props *restrict stars){
for (i = 0; i < stars->lifetimes.n_z; i++) {
for (j = 0; j < stars->lifetimes.n_mass; j++) {
- //index = row_major_index_2d(i,j,stars->lifetimes.n_z,stars->lifetimes.n_mass);
- //stars->lifetimes.dyingtime[index] = log10(temptime[i][j]);
+ // index =
+ // row_major_index_2d(i,j,stars->lifetimes.n_z,stars->lifetimes.n_mass);
+ // stars->lifetimes.dyingtime[index] = log10(temptime[i][j]);
stars->lifetimes.dyingtime[i][j] = log10(temptime[i][j]);
}
}
H5Fclose(file_id);
-
#endif
}
-inline static void allocate_yield_tables(struct stars_props *restrict stars){
+inline static void allocate_yield_tables(struct stars_props *restrict stars) {
/* Allocate SNIa arrays */
- if (posix_memalign((void **)&stars->yields_SNIa, SWIFT_STRUCT_ALIGNMENT, stars->SNIa_n_elements * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yields_SNIa, SWIFT_STRUCT_ALIGNMENT,
+ stars->SNIa_n_elements * sizeof(double)) != 0) {
error("Failed to allocate SNIa yields array");
}
- if (posix_memalign((void **)&stars->yield_SNIa_SPH, SWIFT_STRUCT_ALIGNMENT, stars->SNIa_n_elements * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_SNIa_SPH, SWIFT_STRUCT_ALIGNMENT,
+ stars->SNIa_n_elements * sizeof(double)) != 0) {
error("Failed to allocate SNIa SPH yields array");
}
/* Allocate AGB arrays */
- if (posix_memalign((void **)&stars->yield_AGB.mass, SWIFT_STRUCT_ALIGNMENT, stars->AGB_n_mass * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_AGB.mass, SWIFT_STRUCT_ALIGNMENT,
+ stars->AGB_n_mass * sizeof(double)) != 0) {
error("Failed to allocate AGB mass array");
}
- if (posix_memalign((void **)&stars->yield_AGB.metallicity, SWIFT_STRUCT_ALIGNMENT, stars->AGB_n_z * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_AGB.metallicity,
+ SWIFT_STRUCT_ALIGNMENT,
+ stars->AGB_n_z * sizeof(double)) != 0) {
error("Failed to allocate AGB metallicity array");
}
- if (posix_memalign((void **)&stars->yield_AGB.SPH, SWIFT_STRUCT_ALIGNMENT, stars->AGB_n_z * n_mass_bins * stars->AGB_n_elements * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_AGB.SPH, SWIFT_STRUCT_ALIGNMENT,
+ stars->AGB_n_z * n_mass_bins * stars->AGB_n_elements *
+ sizeof(double)) != 0) {
error("Failed to allocate AGB SPH array");
}
- if (posix_memalign((void **)&stars->yield_AGB.yield, SWIFT_STRUCT_ALIGNMENT, stars->AGB_n_z * stars->AGB_n_mass * stars->AGB_n_elements * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_AGB.yield, SWIFT_STRUCT_ALIGNMENT,
+ stars->AGB_n_z * stars->AGB_n_mass *
+ stars->AGB_n_elements * sizeof(double)) != 0) {
error("Failed to allocate AGB yield array");
}
- if (posix_memalign((void **)&stars->yield_AGB.ejecta_SPH, SWIFT_STRUCT_ALIGNMENT, stars->AGB_n_z * n_mass_bins * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_AGB.ejecta_SPH,
+ SWIFT_STRUCT_ALIGNMENT,
+ stars->AGB_n_z * n_mass_bins * sizeof(double)) != 0) {
error("Failed to allocate AGB ejecta SPH array");
}
- if (posix_memalign((void **)&stars->yield_AGB.ejecta, SWIFT_STRUCT_ALIGNMENT, stars->AGB_n_z * stars->AGB_n_mass * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_AGB.ejecta, SWIFT_STRUCT_ALIGNMENT,
+ stars->AGB_n_z * stars->AGB_n_mass * sizeof(double)) !=
+ 0) {
error("Failed to allocate AGB ejecta array");
}
- if (posix_memalign((void **)&stars->yield_AGB.total_metals_SPH, SWIFT_STRUCT_ALIGNMENT, stars->AGB_n_z * n_mass_bins * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_AGB.total_metals_SPH,
+ SWIFT_STRUCT_ALIGNMENT,
+ stars->AGB_n_z * n_mass_bins * sizeof(double)) != 0) {
error("Failed to allocate AGB total metals SPH array");
}
- if (posix_memalign((void **)&stars->yield_AGB.total_metals, SWIFT_STRUCT_ALIGNMENT, stars->AGB_n_z * stars->AGB_n_mass * sizeof(double)) !=0) {
+ if (posix_memalign(
+ (void **)&stars->yield_AGB.total_metals, SWIFT_STRUCT_ALIGNMENT,
+ stars->AGB_n_z * stars->AGB_n_mass * sizeof(double)) != 0) {
error("Failed to allocate AGB total metals array");
}
/* Allocate SNII arrays */
- if (posix_memalign((void **)&stars->yield_SNII.mass, SWIFT_STRUCT_ALIGNMENT, stars->SNII_n_mass * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_SNII.mass, SWIFT_STRUCT_ALIGNMENT,
+ stars->SNII_n_mass * sizeof(double)) != 0) {
error("Failed to allocate SNII mass array");
}
- if (posix_memalign((void **)&stars->yield_SNII.metallicity, SWIFT_STRUCT_ALIGNMENT, stars->SNII_n_z * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_SNII.metallicity,
+ SWIFT_STRUCT_ALIGNMENT,
+ stars->SNII_n_z * sizeof(double)) != 0) {
error("Failed to allocate SNII metallicity array");
}
- if (posix_memalign((void **)&stars->yield_SNII.SPH, SWIFT_STRUCT_ALIGNMENT, stars->SNII_n_z * n_mass_bins * stars->SNII_n_elements * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_SNII.SPH, SWIFT_STRUCT_ALIGNMENT,
+ stars->SNII_n_z * n_mass_bins * stars->SNII_n_elements *
+ sizeof(double)) != 0) {
error("Failed to allocate SNII SPH array");
}
- if (posix_memalign((void **)&stars->yield_SNII.yield, SWIFT_STRUCT_ALIGNMENT, stars->SNII_n_z * stars->SNII_n_mass * stars->SNII_n_elements * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_SNII.yield, SWIFT_STRUCT_ALIGNMENT,
+ stars->SNII_n_z * stars->SNII_n_mass *
+ stars->SNII_n_elements * sizeof(double)) != 0) {
error("Failed to allocate SNII yield array");
}
- if (posix_memalign((void **)&stars->yield_SNII.ejecta_SPH, SWIFT_STRUCT_ALIGNMENT, stars->SNII_n_z * n_mass_bins * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_SNII.ejecta_SPH,
+ SWIFT_STRUCT_ALIGNMENT,
+ stars->SNII_n_z * n_mass_bins * sizeof(double)) != 0) {
error("Failed to allocate SNII ejecta SPH array");
}
- if (posix_memalign((void **)&stars->yield_SNII.ejecta, SWIFT_STRUCT_ALIGNMENT, stars->SNII_n_z * stars->SNII_n_mass * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_SNII.ejecta, SWIFT_STRUCT_ALIGNMENT,
+ stars->SNII_n_z * stars->SNII_n_mass * sizeof(double)) !=
+ 0) {
error("Failed to allocate SNII ejecta array");
}
- if (posix_memalign((void **)&stars->yield_SNII.total_metals_SPH, SWIFT_STRUCT_ALIGNMENT, stars->SNII_n_z * n_mass_bins * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_SNII.total_metals_SPH,
+ SWIFT_STRUCT_ALIGNMENT,
+ stars->SNII_n_z * n_mass_bins * sizeof(double)) != 0) {
error("Failed to allocate SNII total metals SPH array");
}
- if (posix_memalign((void **)&stars->yield_SNII.total_metals, SWIFT_STRUCT_ALIGNMENT, stars->SNII_n_z * stars->SNII_n_mass * sizeof(double)) !=0) {
+ if (posix_memalign(
+ (void **)&stars->yield_SNII.total_metals, SWIFT_STRUCT_ALIGNMENT,
+ stars->SNII_n_z * stars->SNII_n_mass * sizeof(double)) != 0) {
error("Failed to allocate SNII total metals array");
}
/* Allocate Lifetime arrays */
- if (posix_memalign((void **)&stars->lifetimes.mass, SWIFT_STRUCT_ALIGNMENT, stars->lifetimes.n_mass * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->lifetimes.mass, SWIFT_STRUCT_ALIGNMENT,
+ stars->lifetimes.n_mass * sizeof(double)) != 0) {
error("Failed to allocate lifetime mass array");
}
- if (posix_memalign((void **)&stars->lifetimes.metallicity, SWIFT_STRUCT_ALIGNMENT, stars->lifetimes.n_z * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->lifetimes.metallicity,
+ SWIFT_STRUCT_ALIGNMENT,
+ stars->lifetimes.n_z * sizeof(double)) != 0) {
error("Failed to allocate lifetime metallicity array");
}
- stars->lifetimes.dyingtime = (double **)malloc(stars->lifetimes.n_z * sizeof(double *));
- for(int i = 0; i < stars->lifetimes.n_z; i++){
- stars->lifetimes.dyingtime[i] = (double *)malloc(stars->lifetimes.n_mass * sizeof(double));
+ stars->lifetimes.dyingtime =
+ (double **)malloc(stars->lifetimes.n_z * sizeof(double *));
+ for (int i = 0; i < stars->lifetimes.n_z; i++) {
+ stars->lifetimes.dyingtime[i] =
+ (double *)malloc(stars->lifetimes.n_mass * sizeof(double));
}
/* Allocate SNII factor array */
- if (posix_memalign((void **)&stars->typeII_factor, SWIFT_STRUCT_ALIGNMENT, chemistry_element_count * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->typeII_factor, SWIFT_STRUCT_ALIGNMENT,
+ chemistry_element_count * sizeof(double)) != 0) {
error("Failed to allocate SNII factor array");
}
/* Allocate element name arrays */
- stars->SNIa_element_names = (char **)malloc(stars->SNIa_n_elements * sizeof(char*));
- for(int i = 0; i < stars->SNIa_n_elements; i++){
- stars->SNIa_element_names[i] = (char *)malloc(stars->element_name_length * sizeof(char));
+ stars->SNIa_element_names =
+ (char **)malloc(stars->SNIa_n_elements * sizeof(char *));
+ for (int i = 0; i < stars->SNIa_n_elements; i++) {
+ stars->SNIa_element_names[i] =
+ (char *)malloc(stars->element_name_length * sizeof(char));
}
- stars->SNII_element_names = (char **)malloc(stars->SNII_n_elements * sizeof(char*));
- for(int i = 0; i < stars->SNII_n_elements; i++){
- stars->SNII_element_names[i] = (char *)malloc(stars->element_name_length * sizeof(char));
+ stars->SNII_element_names =
+ (char **)malloc(stars->SNII_n_elements * sizeof(char *));
+ for (int i = 0; i < stars->SNII_n_elements; i++) {
+ stars->SNII_element_names[i] =
+ (char *)malloc(stars->element_name_length * sizeof(char));
}
- stars->AGB_element_names = (char **)malloc(stars->AGB_n_elements * sizeof(char*));
- for(int i = 0; i < stars->AGB_n_elements; i++){
- stars->AGB_element_names[i] = (char *)malloc(stars->element_name_length * sizeof(char));
+ stars->AGB_element_names =
+ (char **)malloc(stars->AGB_n_elements * sizeof(char *));
+ for (int i = 0; i < stars->AGB_n_elements; i++) {
+ stars->AGB_element_names[i] =
+ (char *)malloc(stars->element_name_length * sizeof(char));
}
-
+
/* Allocate array of mass bins for yield and ejecta calculation */
- if (posix_memalign((void **)&stars->yield_mass_bins, SWIFT_STRUCT_ALIGNMENT, n_mass_bins * sizeof(double)) !=0) {
+ if (posix_memalign((void **)&stars->yield_mass_bins, SWIFT_STRUCT_ALIGNMENT,
+ n_mass_bins * sizeof(double)) != 0) {
error("Failed to allocate yield mass bins array");
}
}
-inline static double interpolate_1D(double *array_x, double *array_y, int size, double value){
-
+inline static double interpolate_1D(double *array_x, double *array_y, int size,
+ double value) {
+
double result;
- if (value < array_x[0])
- error("interpolating value less than array min. value %.5e array min %.5e", value, array_x[0]);
- else if (value > array_x[size-1])
- error("interpolating value greater than array max. value %.5e array max %.5e", value, array_x[size-1]);
+ if (value < array_x[0])
+ error("interpolating value less than array min. value %.5e array min %.5e",
+ value, array_x[0]);
+ else if (value > array_x[size - 1])
+ error(
+ "interpolating value greater than array max. value %.5e array max %.5e",
+ value, array_x[size - 1]);
else {
int index = 0;
while (array_x[index] <= value) index++;
- double offset = (array_x[index] - value)/(array_x[index] - array_x[index-1]);
- result = offset*array_y[index-1] + (1 - offset)*array_y[index];
+ double offset =
+ (array_x[index] - value) / (array_x[index] - array_x[index - 1]);
+ result = offset * array_y[index - 1] + (1 - offset) * array_y[index];
}
return result;
}
-inline static void compute_yields(struct stars_props *restrict stars){
+inline static void compute_yields(struct stars_props *restrict stars) {
int index, index_2d;
gsl_interp_accel *accel_ptr;
@@ -480,7 +537,8 @@ inline static void compute_yields(struct stars_props *restrict stars){
}
for (int i = 0; i < stars->SNII_n_z; i++) {
if (stars->yield_SNII.metallicity[i] > 0) {
- stars->yield_SNII.metallicity[i] = log10(stars->yield_SNII.metallicity[i]);
+ stars->yield_SNII.metallicity[i] =
+ log10(stars->yield_SNII.metallicity[i]);
} else {
stars->yield_SNII.metallicity[i] = log_min_metallicity;
}
@@ -508,46 +566,65 @@ inline static void compute_yields(struct stars_props *restrict stars){
/* Loop over elements tracked in EAGLE */
int element_index = 0;
// ALEXEI: better name for eagle_elem?
- for (enum chemistry_element eagle_elem = chemistry_element_H; eagle_elem < chemistry_element_count; eagle_elem++) {
+ for (enum chemistry_element eagle_elem = chemistry_element_H;
+ eagle_elem < chemistry_element_count; eagle_elem++) {
/* SNIa */
- element_index = get_element_index(chemistry_get_element_name(eagle_elem), stars->SNIa_element_names, stars->SNIa_n_elements);
+ element_index =
+ get_element_index(chemistry_get_element_name(eagle_elem),
+ stars->SNIa_element_names, stars->SNIa_n_elements);
stars->yield_SNIa_SPH[eagle_elem] = stars->yields_SNIa[element_index];
/* SNII */
- element_index = get_element_index(chemistry_get_element_name(eagle_elem), stars->SNII_element_names, stars->SNII_n_elements);
+ element_index =
+ get_element_index(chemistry_get_element_name(eagle_elem),
+ stars->SNII_element_names, stars->SNII_n_elements);
for (int i = 0; i < stars->SNII_n_z; i++) {
for (int j = 0; j < stars->SNII_n_mass; j++) {
- index = row_major_index_3d(i, element_index, j, stars->SNII_n_z, stars->SNII_n_elements, stars->SNII_n_mass);
- SNII_yield[j] = stars->yield_SNII.yield[index] * exp(M_LN10 * (-stars->yield_SNII.mass[j]));
+ index = row_major_index_3d(i, element_index, j, stars->SNII_n_z,
+ stars->SNII_n_elements, stars->SNII_n_mass);
+ SNII_yield[j] = stars->yield_SNII.yield[index] *
+ exp(M_LN10 * (-stars->yield_SNII.mass[j]));
}
- gsl_spline_init(SNII_spline_ptr, stars->yield_SNII.mass, SNII_yield, stars->SNII_n_mass);
+ gsl_spline_init(SNII_spline_ptr, stars->yield_SNII.mass, SNII_yield,
+ stars->SNII_n_mass);
for (int k = 0; k < n_mass_bins; k++) {
if (stars->yield_mass_bins[k] < stars->yield_SNII.mass[0])
result = SNII_yield[0];
- else if (stars->yield_mass_bins[k] > stars->yield_SNII.mass[stars->SNII_n_mass - 1])
+ else if (stars->yield_mass_bins[k] >
+ stars->yield_SNII.mass[stars->SNII_n_mass - 1])
result = SNII_yield[stars->SNII_n_mass - 1];
else {
- // Not working with gsl. Aborts inside the evaluation function supposedly because we're trying to interpolate a value out of bounds. Works if SNIIs changed to AGBs though. use own interpolation function for now.
- //result =
- // gsl_spline_eval(SNII_spline_ptr, stars->yield_mass_bins[k], accel_ptr);
- result = interpolate_1D(stars->yield_SNII.mass,SNII_yield, stars->SNII_n_mass,stars->yield_mass_bins[k]);
- }
-
- index = row_major_index_3d(i,eagle_elem,k,stars->SNII_n_z,chemistry_element_count,n_mass_bins);
- stars->yield_SNII.SPH[index] = exp(M_LN10 * stars->yield_mass_bins[k]) * result;
+ // Not working with gsl. Aborts inside the evaluation function
+ // supposedly because we're trying to interpolate a value out of
+ // bounds. Works if SNIIs changed to AGBs though. use own
+ // interpolation function for now.
+ // result =
+ // gsl_spline_eval(SNII_spline_ptr, stars->yield_mass_bins[k],
+ // accel_ptr);
+ result =
+ interpolate_1D(stars->yield_SNII.mass, SNII_yield,
+ stars->SNII_n_mass, stars->yield_mass_bins[k]);
+ }
+
+ index = row_major_index_3d(i, eagle_elem, k, stars->SNII_n_z,
+ chemistry_element_count, n_mass_bins);
+ stars->yield_SNII.SPH[index] =
+ exp(M_LN10 * stars->yield_mass_bins[k]) * result;
}
}
for (int i = 0; i < stars->SNII_n_z; i++) {
for (int k = 0; k < n_mass_bins; k++) {
- index_2d = row_major_index_2d(i,k,stars->SNII_n_z,n_mass_bins);
- index = row_major_index_3d(i,eagle_elem,k,stars->SNII_n_z,chemistry_element_count,n_mass_bins);
+ index_2d = row_major_index_2d(i, k, stars->SNII_n_z, n_mass_bins);
+ index = row_major_index_3d(i, eagle_elem, k, stars->SNII_n_z,
+ chemistry_element_count, n_mass_bins);
if (strcmp(chemistry_get_element_name(eagle_elem), "Hydrogen") != 0 ||
strcmp(chemistry_get_element_name(eagle_elem), "Helium") != 0) {
stars->yield_SNII.total_metals_SPH[index_2d] +=
- (stars->typeII_factor[eagle_elem] - 1) * stars->yield_SNII.SPH[index];
+ (stars->typeII_factor[eagle_elem] - 1) *
+ stars->yield_SNII.SPH[index];
}
stars->yield_SNII.SPH[index] *= stars->typeII_factor[eagle_elem];
@@ -555,35 +632,47 @@ inline static void compute_yields(struct stars_props *restrict stars){
}
/* AGB */
- element_index = get_element_index(chemistry_get_element_name(eagle_elem), stars->AGB_element_names, stars->AGB_n_elements);
+ element_index =
+ get_element_index(chemistry_get_element_name(eagle_elem),
+ stars->AGB_element_names, stars->AGB_n_elements);
if (element_index < 0) {
for (int i = 0; i < stars->AGB_n_z; i++) {
for (int j = 0; j < n_mass_bins; j++) {
- index = row_major_index_3d(i,eagle_elem,j,stars->AGB_n_z,chemistry_element_count,n_mass_bins);
+ index = row_major_index_3d(i, eagle_elem, j, stars->AGB_n_z,
+ chemistry_element_count, n_mass_bins);
stars->yield_AGB.SPH[index] = 0.0;
- }
+ }
}
} else {
for (int i = 0; i < stars->AGB_n_z; i++) {
for (int j = 0; j < stars->AGB_n_mass; j++) {
- index = row_major_index_3d(i, element_index, j, stars->AGB_n_z, stars->AGB_n_elements, stars->AGB_n_mass);
- AGB_yield[j] = stars->yield_AGB.yield[index] * exp(M_LN10 * (-stars->yield_AGB.mass[j]));
+ index = row_major_index_3d(i, element_index, j, stars->AGB_n_z,
+ stars->AGB_n_elements, stars->AGB_n_mass);
+ AGB_yield[j] = stars->yield_AGB.yield[index] *
+ exp(M_LN10 * (-stars->yield_AGB.mass[j]));
}
- gsl_spline_init(AGB_spline_ptr, stars->yield_AGB.mass, AGB_yield, stars->AGB_n_mass);
+ gsl_spline_init(AGB_spline_ptr, stars->yield_AGB.mass, AGB_yield,
+ stars->AGB_n_mass);
for (int j = 0; j < n_mass_bins; j++) {
if (stars->yield_mass_bins[j] < stars->yield_AGB.mass[0])
result = AGB_yield[0];
- else if (stars->yield_mass_bins[j] > stars->yield_AGB.mass[stars->AGB_n_mass - 1])
+ else if (stars->yield_mass_bins[j] >
+ stars->yield_AGB.mass[stars->AGB_n_mass - 1])
result = AGB_yield[stars->AGB_n_mass - 1];
- else
- result =
- gsl_spline_eval(AGB_spline_ptr, stars->yield_mass_bins[j], accel_ptr);
-
- index = row_major_index_3d(i,eagle_elem,j,stars->AGB_n_z,chemistry_element_count,n_mass_bins);
- stars->yield_AGB.SPH[index] = exp(M_LN10 * stars->yield_mass_bins[j]) * result;
- if (element_index == 10) message("j %d AGB SPH %.5e result %.5e exponential %.5e", j, stars->yield_AGB.SPH[index], result, exp(M_LN10 * stars->yield_mass_bins[j]) );
+ else
+ result = gsl_spline_eval(AGB_spline_ptr, stars->yield_mass_bins[j],
+ accel_ptr);
+
+ index = row_major_index_3d(i, eagle_elem, j, stars->AGB_n_z,
+ chemistry_element_count, n_mass_bins);
+ stars->yield_AGB.SPH[index] =
+ exp(M_LN10 * stars->yield_mass_bins[j]) * result;
+ if (element_index == 10)
+ message("j %d AGB SPH %.5e result %.5e exponential %.5e", j,
+ stars->yield_AGB.SPH[index], result,
+ exp(M_LN10 * stars->yield_mass_bins[j]));
}
}
}
@@ -591,68 +680,76 @@ inline static void compute_yields(struct stars_props *restrict stars){
gsl_spline_free(SNII_spline_ptr);
gsl_spline_free(AGB_spline_ptr);
gsl_interp_accel_free(accel_ptr);
-
}
inline static void compute_ejecta(struct stars_props *restrict stars) {
-
+
gsl_interp_accel *accel_ptr;
gsl_spline *SNII_spline_ptr, *AGB_spline_ptr;
- // Do we really need SNII_yield and AGB_yield, they're not used simultaneously, so can use only one?
+ // Do we really need SNII_yield and AGB_yield, they're not used
+ // simultaneously, so can use only one?
double SNII_yield[stars->SNII_n_mass];
double AGB_yield[stars->AGB_n_mass];
float result;
-
- accel_ptr = gsl_interp_accel_alloc();
- SNII_spline_ptr = gsl_spline_alloc(gsl_interp_linear, stars->SNII_n_mass);
+
+ accel_ptr = gsl_interp_accel_alloc();
+ SNII_spline_ptr = gsl_spline_alloc(gsl_interp_linear, stars->SNII_n_mass);
int index;
-
+
for (int i = 0; i < stars->SNII_n_z; i++) {
for (int k = 0; k < stars->SNII_n_mass; k++) {
- index = row_major_index_2d(i,k,stars->SNII_n_z,stars->SNII_n_mass);
- SNII_yield[k] = stars->yield_SNII.ejecta[index] * exp(M_LN10 * (-stars->yield_SNII.mass[k]));
+ index = row_major_index_2d(i, k, stars->SNII_n_z, stars->SNII_n_mass);
+ SNII_yield[k] = stars->yield_SNII.ejecta[index] *
+ exp(M_LN10 * (-stars->yield_SNII.mass[k]));
}
- gsl_spline_init(SNII_spline_ptr, stars->yield_SNII.mass, SNII_yield, stars->SNII_n_mass);
-
- for (int k = 0; k < n_mass_bins; k++) {
- if (stars->yield_mass_bins[k] < stars->yield_SNII.mass[0])
+ gsl_spline_init(SNII_spline_ptr, stars->yield_SNII.mass, SNII_yield,
+ stars->SNII_n_mass);
+
+ for (int k = 0; k < n_mass_bins; k++) {
+ if (stars->yield_mass_bins[k] < stars->yield_SNII.mass[0])
result = SNII_yield[0];
- else if (stars->yield_mass_bins[k] > stars->yield_SNII.mass[stars->SNII_n_mass - 1])
+ else if (stars->yield_mass_bins[k] >
+ stars->yield_SNII.mass[stars->SNII_n_mass - 1])
result = SNII_yield[stars->SNII_n_mass - 1];
else
- result =
- gsl_spline_eval(SNII_spline_ptr, stars->yield_mass_bins[k], accel_ptr);
-
- index = row_major_index_2d(i,k,stars->SNII_n_z,stars->SNII_n_mass);
- stars->yield_SNII.ejecta_SPH[index] = exp(M_LN10 * stars->yield_mass_bins[k]) * result;
- }
- }
-
- for (int i = 0; i < stars->SNII_n_z; i++) {
+ result = gsl_spline_eval(SNII_spline_ptr, stars->yield_mass_bins[k],
+ accel_ptr);
+
+ index = row_major_index_2d(i, k, stars->SNII_n_z, stars->SNII_n_mass);
+ stars->yield_SNII.ejecta_SPH[index] =
+ exp(M_LN10 * stars->yield_mass_bins[k]) * result;
+ }
+ }
+
+ for (int i = 0; i < stars->SNII_n_z; i++) {
for (int k = 0; k < stars->SNII_n_mass; k++) {
- index = row_major_index_2d(i,k,stars->SNII_n_z,stars->SNII_n_mass);
- SNII_yield[k] = stars->yield_SNII.total_metals[index] * exp(M_LN10 * (-stars->yield_SNII.mass[k]));
+ index = row_major_index_2d(i, k, stars->SNII_n_z, stars->SNII_n_mass);
+ SNII_yield[k] = stars->yield_SNII.total_metals[index] *
+ exp(M_LN10 * (-stars->yield_SNII.mass[k]));
}
-
- gsl_spline_init(SNII_spline_ptr, stars->yield_SNII.mass, SNII_yield, stars->SNII_n_mass);
-
+
+ gsl_spline_init(SNII_spline_ptr, stars->yield_SNII.mass, SNII_yield,
+ stars->SNII_n_mass);
+
for (int k = 0; k < n_mass_bins; k++) {
- if (stars->yield_mass_bins[k] < stars->yield_SNII.mass[0])
- result = SNII_yield[0];
- else if (stars->yield_mass_bins[k] > stars->yield_SNII.mass[stars->SNII_n_mass - 1])
- result = SNII_yield[stars->SNII_n_mass - 1];
- else
- result =
- gsl_spline_eval(SNII_spline_ptr, stars->yield_mass_bins[k], accel_ptr);
-
- index = row_major_index_2d(i,k,stars->SNII_n_z,stars->SNII_n_mass);
- stars->yield_SNII.total_metals_SPH[index] = exp(M_LN10 * stars->yield_mass_bins[k]) * result;
- }
- }
+ if (stars->yield_mass_bins[k] < stars->yield_SNII.mass[0])
+ result = SNII_yield[0];
+ else if (stars->yield_mass_bins[k] >
+ stars->yield_SNII.mass[stars->SNII_n_mass - 1])
+ result = SNII_yield[stars->SNII_n_mass - 1];
+ else
+ result = gsl_spline_eval(SNII_spline_ptr, stars->yield_mass_bins[k],
+ accel_ptr);
+
+ index = row_major_index_2d(i, k, stars->SNII_n_z, stars->SNII_n_mass);
+ stars->yield_SNII.total_metals_SPH[index] =
+ exp(M_LN10 * stars->yield_mass_bins[k]) * result;
+ }
+ }
gsl_spline_free(SNII_spline_ptr);
@@ -661,45 +758,53 @@ inline static void compute_ejecta(struct stars_props *restrict stars) {
for (int i = 0; i < stars->AGB_n_z; i++) {
for (int k = 0; k < stars->AGB_n_mass; k++) {
- index = row_major_index_2d(i,k,stars->SNII_n_z,stars->SNII_n_mass);
- AGB_yield[k] = stars->yield_AGB.ejecta[index] / exp(M_LN10 * stars->yield_AGB.mass[k]);
+ index = row_major_index_2d(i, k, stars->SNII_n_z, stars->SNII_n_mass);
+ AGB_yield[k] = stars->yield_AGB.ejecta[index] /
+ exp(M_LN10 * stars->yield_AGB.mass[k]);
}
- gsl_spline_init(AGB_spline_ptr, stars->yield_AGB.mass, AGB_yield, stars->AGB_n_mass);
+ gsl_spline_init(AGB_spline_ptr, stars->yield_AGB.mass, AGB_yield,
+ stars->AGB_n_mass);
for (int k = 0; k < n_mass_bins; k++) {
if (stars->yield_mass_bins[k] < stars->yield_AGB.mass[0])
result = AGB_yield[0];
- else if (stars->yield_mass_bins[k] > stars->yield_AGB.mass[stars->AGB_n_mass - 1])
+ else if (stars->yield_mass_bins[k] >
+ stars->yield_AGB.mass[stars->AGB_n_mass - 1])
result = AGB_yield[stars->AGB_n_mass - 1];
else
- result =
- gsl_spline_eval(AGB_spline_ptr, stars->yield_mass_bins[k], accel_ptr);
+ result = gsl_spline_eval(AGB_spline_ptr, stars->yield_mass_bins[k],
+ accel_ptr);
- index = row_major_index_2d(i,k,stars->SNII_n_z,stars->SNII_n_mass);
- stars->yield_AGB.ejecta_SPH[index] = exp(M_LN10 * stars->yield_mass_bins[k]) * result;
+ index = row_major_index_2d(i, k, stars->SNII_n_z, stars->SNII_n_mass);
+ stars->yield_AGB.ejecta_SPH[index] =
+ exp(M_LN10 * stars->yield_mass_bins[k]) * result;
}
}
for (int i = 0; i < stars->AGB_n_z; i++) {
for (int k = 0; k < stars->AGB_n_mass; k++) {
- index = row_major_index_2d(i,k,stars->SNII_n_z,stars->SNII_n_mass);
- AGB_yield[k] = stars->yield_AGB.total_metals[index] * exp(M_LN10 * (-stars->yield_AGB.mass[k]));
+ index = row_major_index_2d(i, k, stars->SNII_n_z, stars->SNII_n_mass);
+ AGB_yield[k] = stars->yield_AGB.total_metals[index] *
+ exp(M_LN10 * (-stars->yield_AGB.mass[k]));
}
- gsl_spline_init(AGB_spline_ptr, stars->yield_AGB.mass, AGB_yield, stars->AGB_n_mass);
+ gsl_spline_init(AGB_spline_ptr, stars->yield_AGB.mass, AGB_yield,
+ stars->AGB_n_mass);
for (int k = 0; k < n_mass_bins; k++) {
if (stars->yield_mass_bins[k] < stars->yield_AGB.mass[0])
result = AGB_yield[0];
- else if (stars->yield_mass_bins[k] > stars->yield_AGB.mass[stars->AGB_n_mass - 1])
+ else if (stars->yield_mass_bins[k] >
+ stars->yield_AGB.mass[stars->AGB_n_mass - 1])
result = AGB_yield[stars->AGB_n_mass - 1];
else
- result =
- gsl_spline_eval(AGB_spline_ptr, stars->yield_mass_bins[k], accel_ptr);
+ result = gsl_spline_eval(AGB_spline_ptr, stars->yield_mass_bins[k],
+ accel_ptr);
- index = row_major_index_2d(i,k,stars->SNII_n_z,stars->SNII_n_mass);
- stars->yield_AGB.total_metals_SPH[index] = exp(M_LN10 * stars->yield_mass_bins[k]) * result;
+ index = row_major_index_2d(i, k, stars->SNII_n_z, stars->SNII_n_mass);
+ stars->yield_AGB.total_metals_SPH[index] =
+ exp(M_LN10 * stars->yield_mass_bins[k]) * result;
}
}
@@ -707,5 +812,4 @@ inline static void compute_ejecta(struct stars_props *restrict stars) {
gsl_interp_accel_free(accel_ptr);
}
-
#endif
diff --git a/src/stars/const/stars.h b/src/stars/const/stars.h
index aec793e0053fa64d768516629f206ef777e175ab..28debcdaf4852b76c96eaad0c08f8a09cfbcc27f 100644
--- a/src/stars/const/stars.h
+++ b/src/stars/const/stars.h
@@ -75,7 +75,7 @@ __attribute__((always_inline)) INLINE static void stars_init_spart(
* @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) {
+ struct spart* restrict sp, float dt_drift) {
const float h_inv = 1.f / sp->h;
@@ -85,7 +85,6 @@ __attribute__((always_inline)) INLINE static void stars_predict_extra(
sp->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
else
sp->h *= expf(w1);
-
}
/**
@@ -184,32 +183,42 @@ __attribute__((always_inline)) INLINE static void stars_reset_acceleration(
*/
__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) {
-
+ 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;
+ 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;
+ 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.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.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;
+ 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 */
+ /* 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);
+ FILE* feedback_output = fopen("feedback_properties.dat", "w");
+ fprintf(feedback_output, "%.5e \n%.5e \n", total_sn,
+ stars_properties->feedback_timescale);
fclose(feedback_output);
}
@@ -218,11 +227,10 @@ __attribute__((always_inline)) INLINE static void stars_evolve_spart(
/* 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){}
-
+inline static void stars_evolve_init(struct swift_params* params,
+ struct stars_props* restrict stars) {}
/**
* @brief Reset acceleration fields of a particle
@@ -245,5 +253,4 @@ __attribute__((always_inline)) INLINE static void stars_reset_feedback(
#endif
}
-
#endif /* SWIFT_CONST_STARS_H */
diff --git a/src/stars/const/stars_iact.h b/src/stars/const/stars_iact.h
index 4b62e4d5a90e04e01e2848edcd7ea4fde48b7561..c0a2fb3f429043fe3134e6c671f5d5594a63744e 100644
--- a/src/stars/const/stars_iact.h
+++ b/src/stars/const/stars_iact.h
@@ -15,12 +15,11 @@
* @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) {
+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;
@@ -32,7 +31,7 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
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;
@@ -45,8 +44,9 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
/* 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);
+ /* 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 */
@@ -57,7 +57,7 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
}
/**
- * @brief Increases thermal energy of particle due
+ * @brief Increases thermal energy of particle due
* to feedback by specified amount
*
* @param du change in internal energy
@@ -65,15 +65,16 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
* @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) {
+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.
+ // 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;
}
@@ -90,15 +91,15 @@ static inline void thermal_feedback(float du, struct part * restrict p,
* @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
+ * @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) {
+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. */
@@ -111,66 +112,97 @@ runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
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);
+ 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;
+ 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);
+ 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);
- 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;
+ 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;
-
+ 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;
+ 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;
+ 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;
+ 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;
+ 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;
+ 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;
+ 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;
+ 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;
+ 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++) {
@@ -180,37 +212,51 @@ runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
/* 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);
+ 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);
+ // 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;
+ 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;
+ 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);
+ 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));
+ 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));
+ 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
index 62b9caeab54ecbfbe85858a31bcddecd07fd6225..d3f184f9aac6592f5b40b4b9fc917020f9eac6e3 100644
--- a/src/stars/const/stars_io.h
+++ b/src/stars/const/stars_io.h
@@ -97,7 +97,7 @@ INLINE static void stars_props_init(struct stars_props *sp,
const struct unit_system *us,
struct swift_params *params,
const struct hydro_props *p,
- const struct cosmology *cosmo) {
+ const struct cosmology *cosmo) {
/* Kernel properties */
sp->eta_neighbours = parser_get_opt_param_float(
@@ -129,39 +129,53 @@ INLINE static void stars_props_init(struct stars_props *sp,
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);
+ /* 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);
+ 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);
+ 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);
+ 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);
+ 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
+ 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) */
+ /* 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.
+ 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;
-
}
/**
diff --git a/src/stars/const/stars_part.h b/src/stars/const/stars_part.h
index 88a954ee3896ee8029c48ec3c6c71aa926a54d21..5927efff9b3abdb78dcdbf693d4044a53b1c73fa 100644
--- a/src/stars/const/stars_part.h
+++ b/src/stars/const/stars_part.h
@@ -23,7 +23,7 @@
#include <stdlib.h>
/* Read chemistry */
-#include "chemistry_struct.h"
+#include "chemistry_struct.h"
/**
* @brief Particle fields for the star particles.
@@ -83,7 +83,7 @@ struct spart {
/* Mass fractions of ejecta */
struct chemistry_part_data chemistry_data;
-
+
float ejecta_specific_thermal_energy;
/* Number of type 1a SNe per unit mass */
@@ -91,7 +91,8 @@ struct spart {
} to_distribute;
- /* kernel normalisation factor (equivalent to metalweight_norm in eagle_enrich.c:811, TODO: IMPROVE COMMENT) */
+ /* 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 */
@@ -152,7 +153,8 @@ struct stars_props {
/* 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) */
+ /* 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 */
@@ -170,7 +172,8 @@ struct stars_props {
/* 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) */
+ /* 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) */
diff --git a/src/tools.c b/src/tools.c
index 8ce0ff6dd987182852819d48e0ea4e453036f513..ef548c5eff140211abc9b6882222cce074f6b754 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -453,7 +453,8 @@ void pairs_all_stars_density(struct runner *r, struct cell *ci,
/* Hit or miss? */
if (r2 < hig2) {
/* Interact */
- runner_iact_nonsym_stars_density(r2, dx, hi, pj->h, spi, pj, cosmo, stars_properties, xpj, 0);
+ runner_iact_nonsym_stars_density(r2, dx, hi, pj->h, spi, pj, cosmo,
+ stars_properties, xpj, 0);
}
}
}
@@ -484,7 +485,8 @@ void pairs_all_stars_density(struct runner *r, struct cell *ci,
/* Hit or miss? */
if (r2 < hjg2) {
/* Interact */
- runner_iact_nonsym_stars_density(r2, dx, hj, pi->h, spj, pi, cosmo, stars_properties, xpi, 0);
+ runner_iact_nonsym_stars_density(r2, dx, hj, pi->h, spj, pi, cosmo,
+ stars_properties, xpi, 0);
}
}
}
@@ -672,7 +674,8 @@ void self_all_stars_density(struct runner *r, struct cell *ci) {
/* Hit or miss? */
if (r2 > 0.f && r2 < hig2) {
/* Interact */
- runner_iact_nonsym_stars_density(r2, dxi, hi, hj, spi, pj, cosmo, stars_properties, xpj, 0);
+ runner_iact_nonsym_stars_density(r2, dxi, hi, hj, spi, pj, cosmo,
+ stars_properties, xpj, 0);
}
}
}
diff --git a/tests/testFeedback.c b/tests/testFeedback.c
index a09188716aba07829a259259fe5e315801fff0a8..eac57882210610cff739c755b904ba31847c67f2 100644
--- a/tests/testFeedback.c
+++ b/tests/testFeedback.c
@@ -55,7 +55,8 @@ int main(int argc, char *argv[]) {
hydro_props_init(&hydro_properties, &phys_const, &us, params);
/* Init star properties */
- stars_props_init(&stars_properties, &phys_const, &us, params, &hydro_properties, &cosmo);
+ stars_props_init(&stars_properties, &phys_const, &us, params,
+ &hydro_properties, &cosmo);
/* Read yield tables */
stars_evolve_init(params, &stars_properties);
@@ -79,12 +80,16 @@ int main(int argc, char *argv[]) {
"| Ne | Mg | Si | Fe | per solar mass\n");
float Gyr_to_s = 3.154e16;
- float dt = 0.1 * Gyr_to_s / units_cgs_conversion_factor(&us,UNIT_CONV_TIME);
- float max_age = 13.f * Gyr_to_s / units_cgs_conversion_factor(&us,UNIT_CONV_TIME);
+ float dt = 0.1 * Gyr_to_s / units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
+ float max_age =
+ 13.f * Gyr_to_s / units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
for (float age = 0; age <= max_age; age += dt) {
compute_stellar_evolution(&stars_properties, &sp, &us, age, dt);
- float age_Gyr = age * units_cgs_conversion_factor(&us,UNIT_CONV_TIME) / Gyr_to_s;
- fprintf(AGB_output, "%f %e %e ", age_Gyr, sp.to_distribute.mass / sp.mass_init, sp.chemistry_data.metal_mass_fraction_from_AGB / sp.mass_init);
+ float age_Gyr =
+ age * units_cgs_conversion_factor(&us, UNIT_CONV_TIME) / Gyr_to_s;
+ fprintf(AGB_output, "%f %e %e ", age_Gyr,
+ sp.to_distribute.mass / sp.mass_init,
+ sp.chemistry_data.metal_mass_fraction_from_AGB / sp.mass_init);
for (int i = 0; i < chemistry_element_count; i++)
fprintf(AGB_output, "%e ", sp.metals_released[i]);
fprintf(AGB_output, "\n");
diff --git a/tests/testStellarEvolution.c b/tests/testStellarEvolution.c
index 6747287ee3aca7b27bd8e7be91525dc08ad2635e..c9c9e793e76a23fd5102b9bcd998fb1f6417d030 100644
--- a/tests/testStellarEvolution.c
+++ b/tests/testStellarEvolution.c
@@ -21,9 +21,9 @@
#include "../config.h"
#include "hydro.h"
#include "physical_constants.h"
+#include "stars.h"
#include "swift.h"
#include "units.h"
-#include "stars.h"
int main(int argc, char **argv) {
/* Declare relevant structs */
@@ -61,7 +61,8 @@ int main(int argc, char **argv) {
hydro_props_init(&hydro_properties, &phys_const, &us, params);
/* Init star properties */
- stars_props_init(&stars_properties, &phys_const, &us, params, &hydro_properties, &cosmo);
+ stars_props_init(&stars_properties, &phys_const, &us, params,
+ &hydro_properties, &cosmo);
/* Init spart */
stars_first_init_spart(&sp);
@@ -72,18 +73,27 @@ int main(int argc, char **argv) {
stars_evolve_spart(&sp, &stars_properties, &cosmo, &us, current_time, dt);
for (int i = 0; i < 9; i++) {
- message("element %d to distribute fraction %.5e", i, sp.to_distribute.chemistry_data.metal_mass_fraction[i]);
+ message("element %d to distribute fraction %.5e", i,
+ sp.to_distribute.chemistry_data.metal_mass_fraction[i]);
}
- message("to distribute mass %.5e",sp.to_distribute.mass);
- message("to distribute num_SNIa %.5e",sp.to_distribute.num_SNIa);
- message("to distribute metal_mass_fraction_total %.5e",sp.to_distribute.chemistry_data.metal_mass_fraction_total);
- message("to distribute mass_from_AGB %.5e",sp.to_distribute.chemistry_data.mass_from_AGB);
- message("to distribute metal_mass_fraction_from_AGB %.5e",sp.to_distribute.chemistry_data.metal_mass_fraction_from_AGB);
- message("to distribute mass_from_SNII %.5e",sp.to_distribute.chemistry_data.mass_from_SNII);
- message("to distribute metal_mass_fraction_from_SNII %.5e",sp.to_distribute.chemistry_data.metal_mass_fraction_from_SNII);
- message("to distribute mass_from_SNIa %.5e",sp.to_distribute.chemistry_data.mass_from_SNIa);
- message("to distribute metal_mass_fraction_from_SNIa %.5e",sp.to_distribute.chemistry_data.metal_mass_fraction_from_SNIa);
- message("to distribute iron_mass_fraction_from_SNIa %.5e",sp.to_distribute.chemistry_data.iron_mass_fraction_from_SNIa);
+ message("to distribute mass %.5e", sp.to_distribute.mass);
+ message("to distribute num_SNIa %.5e", sp.to_distribute.num_SNIa);
+ message("to distribute metal_mass_fraction_total %.5e",
+ sp.to_distribute.chemistry_data.metal_mass_fraction_total);
+ message("to distribute mass_from_AGB %.5e",
+ sp.to_distribute.chemistry_data.mass_from_AGB);
+ message("to distribute metal_mass_fraction_from_AGB %.5e",
+ sp.to_distribute.chemistry_data.metal_mass_fraction_from_AGB);
+ message("to distribute mass_from_SNII %.5e",
+ sp.to_distribute.chemistry_data.mass_from_SNII);
+ message("to distribute metal_mass_fraction_from_SNII %.5e",
+ sp.to_distribute.chemistry_data.metal_mass_fraction_from_SNII);
+ message("to distribute mass_from_SNIa %.5e",
+ sp.to_distribute.chemistry_data.mass_from_SNIa);
+ message("to distribute metal_mass_fraction_from_SNIa %.5e",
+ sp.to_distribute.chemistry_data.metal_mass_fraction_from_SNIa);
+ message("to distribute iron_mass_fraction_from_SNIa %.5e",
+ sp.to_distribute.chemistry_data.iron_mass_fraction_from_SNIa);
message("done test");