Skip to content
GitLab
Commit 519718d7 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Smoothed metals in enrichment

parent 71570124
Expand all Hide whitespace changes
Inline Side-by-side
examples/SubgridTests/StellarEvolution/plot_box_evolution.py
View file @ 519718d7
......@@ -90,8 +90,14 @@ Msun_in_cgs = 1.98848e33
# Declare arrays to store SWIFT data
swift_box_gas_mass = zeros(n_snapshots)
swift_box_gas_mass_AGB = zeros(n_snapshots)
swift_box_gas_mass_SNII = zeros(n_snapshots)
swift_box_gas_mass_SNIa = zeros(n_snapshots)
swift_box_star_mass = zeros(n_snapshots)
swift_box_gas_metal_mass = zeros(n_snapshots)
swift_box_gas_metal_mass_AGB = zeros(n_snapshots)
swift_box_gas_metal_mass_SNII = zeros(n_snapshots)
swift_box_gas_metal_mass_SNIa = zeros(n_snapshots)
swift_element_mass = zeros((n_snapshots,n_elements))
swift_internal_energy = zeros(n_snapshots)
swift_kinetic_energy = zeros(n_snapshots)
......@@ -108,6 +114,14 @@ for i in range(n_snapshots):
masses = sim["/PartType0/Masses"][:]
swift_box_gas_mass[i] = np.sum(masses)
AGB_mass = sim["/PartType0/MassesFromAGB"][:]
SNII_mass = sim["/PartType0/MassesFromSNII"][:]
SNIa_mass = sim["/PartType0/MassesFromSNIa"][:]
swift_box_gas_mass_AGB[i] = np.sum(AGB_mass)
swift_box_gas_mass_SNII[i] = np.sum(SNII_mass)
swift_box_gas_mass_SNIa[i] = np.sum(SNIa_mass)
Z_star = sim["/PartType4/MetalMassFractions"][0]
star_masses = sim["/PartType4/Masses"][:]
swift_box_star_mass[i] = np.sum(star_masses)
......@@ -115,6 +129,14 @@ for i in range(n_snapshots):
metallicities = sim["/PartType0/MetalMassFractions"][:]
swift_box_gas_metal_mass[i] = np.sum(metallicities * masses)
AGB_Z_frac = sim["/PartType0/MetalMassFractionsFromAGB"][:]
SNII_Z_frac = sim["/PartType0/MetalMassFractionsFromSNII"][:]
SNIa_Z_frac = sim["/PartType0/MetalMassFractionsFromSNIa"][:]
swift_box_gas_metal_mass_AGB[i] = np.sum(AGB_Z_frac * masses)
swift_box_gas_metal_mass_SNII[i] = np.sum(SNII_Z_frac * masses)
swift_box_gas_metal_mass_SNIa[i] = np.sum(SNIa_Z_frac * masses)
element_abundances = sim["/PartType0/ElementMassFractions"][:][:]
for j in range(n_elements):
swift_element_mass[i,j] = np.sum(element_abundances[:,j] * masses)
......@@ -136,58 +158,32 @@ for i in range(n_snapshots):
filename = "./StellarEvolutionSolution/Z_%.4f/StellarEvolutionTotal.txt"%Z_star
# Read EAGLE test output
with open(filename) as f:
eagle_categories = f.readline()
eagle_data = f.readlines()
eagle_data = [x.strip() for x in eagle_data]
# Declare arrays to store EAGLE test output
eagle_time_Gyr = zeros(len(eagle_data))
eagle_total_mass = zeros(len(eagle_data))
eagle_total_metal_mass = zeros(len(eagle_data))
eagle_total_element_mass = zeros((len(eagle_data),n_elements))
eagle_energy_from_mass_cgs = zeros(len(eagle_data))
eagle_energy_ejecta_cgs = zeros(len(eagle_data))
# Populate arrays with data from EAGLE test output
i = 0
for line in eagle_data:
enrich_to_date = line.split(' ')
eagle_time_Gyr[i] = float(enrich_to_date[0])
eagle_total_mass[i] = float(enrich_to_date[1]) * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
eagle_total_metal_mass[i] = float(enrich_to_date[2]) * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
for j in range(n_elements):
eagle_total_element_mass[i,j] = float(enrich_to_date[3+j]) * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
eagle_energy_from_mass_cgs[i] = eagle_total_mass[i] * Msun_in_cgs * swift_mean_u_start * unit_int_energy_in_cgs
eagle_energy_ejecta_cgs[i] = 0.5 * (eagle_total_mass[i] * Msun_in_cgs) * ejecta_vel_cgs**2
i += 1
data = loadtxt(filename)
eagle_time_Gyr = data[:,0]
eagle_total_mass = data[:,1] * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
eagle_total_metal_mass = data[:,2] * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
eagle_total_element_mass = data[:, 3:3+n_elements] * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
eagle_energy_from_mass_cgs = eagle_total_mass * Msun_in_cgs * swift_mean_u_start * unit_int_energy_in_cgs
eagle_energy_ejecta_cgs = 0.5 * (eagle_total_mass * Msun_in_cgs) * ejecta_vel_cgs**2
# Read the mass per channel
filename = "./StellarEvolutionSolution/Z_%.4f/StellarEvolutionAGB.txt"%Z_star
data = loadtxt(filename)
eagle_total_mass_AGB = data[:,1] * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
eagle_total_metal_mass_AGB = data[:,2] * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
filename = "./StellarEvolutionSolution/Z_%.4f/StellarEvolutionII.txt"%Z_star
data = loadtxt(filename)
eagle_total_mass_SNII = data[:,1] * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
eagle_total_metal_mass_SNII = data[:,2] * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
# Read the number of SNIa
filename = "./StellarEvolutionSolution/Z_%.4f/StellarEvolutionIa.txt"%Z_star
with open(filename) as f:
eagle_categories = f.readline()
eagle_data = f.readlines()
i = 0
N_SNIa = zeros(len(eagle_data))
for line in eagle_data:
enrich_to_date = line.split(' ')
N_SNIa[i] = float(enrich_to_date[-2]) * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
i += 1
cumulative_N_SNIa = np.cumsum(N_SNIa)
eagle_energy_SNIa_cgs = cumulative_N_SNIa * E_SNIa_cgs
# SNII energy
N_SNII = 0.017362 * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
eagle_energy_SNII_cgs = np.ones(len(eagle_data)) * N_SNII * E_SNII_cgs
eagle_energy_SNII_cgs[eagle_time_Gyr < 0.03] = 0.
# Total energy
eagle_energy_total_cgs = eagle_energy_ejecta_cgs + eagle_energy_from_mass_cgs + eagle_energy_SNIa_cgs
data = loadtxt(filename)
eagle_total_mass_SNIa = data[:,1] * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
eagle_total_metal_mass_SNIa = data[:,2] * stellar_mass / Msun_in_cgs * unit_mass_in_cgs
# Plot the interesting quantities
figure()
......@@ -195,19 +191,25 @@ suptitle("Star metallicity Z = %.4f"%Z_star)
# Box gas mass --------------------------------
subplot(221)
plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (swift_box_gas_mass[1:] - swift_box_gas_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', marker = "*", ms=0.5, label='swift')
plot(eagle_time_Gyr[1:],eagle_total_mass[:-1],linewidth=0.5,color='r',label='eagle test total', ls='--')
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in total gas particle mass (Msun)", labelpad=2)
plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (swift_box_gas_mass[1:] - swift_box_gas_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', label='Total')
plot(t * unit_time_in_cgs / Gyr_in_cgs, swift_box_gas_mass_AGB * unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='C0', label='AGB')
plot(t * unit_time_in_cgs / Gyr_in_cgs, swift_box_gas_mass_SNII * unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='C1', label='SNII')
plot(t * unit_time_in_cgs / Gyr_in_cgs, swift_box_gas_mass_SNIa * unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='C2', label='SNIa')
plot(eagle_time_Gyr[1:],eagle_total_mass[:-1],linewidth=0.5, color='k', ls='--')
plot(eagle_time_Gyr[1:],eagle_total_mass_AGB[:-1],linewidth=0.5, color='C0', ls='--')
plot(eagle_time_Gyr[1:],eagle_total_mass_SNII[:-1],linewidth=0.5, color='C1', ls='--')
plot(eagle_time_Gyr[1:],eagle_total_mass_SNIa[:-1],linewidth=0.5, color='C2', ls='--')
legend(loc="lower right", ncol=2, fontsize=8)
xlabel("${\\rm Time~[Gyr]}$", labelpad=0)
ylabel("Change in total gas particle mass ${[\\rm M_\\odot]}$", labelpad=2)
ticklabel_format(style='sci', axis='y', scilimits=(0,0))
legend()
# Box star mass --------------------------------
subplot(222)
plot(t * unit_time_in_cgs / Gyr_in_cgs, (swift_box_star_mass)* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', marker = "*", ms=0.5, label='swift')
plot(eagle_time_Gyr[1:], swift_box_star_mass[0] * unit_mass_in_cgs / Msun_in_cgs - eagle_total_mass[:-1],linewidth=0.5,color='r',label='eagle test total')
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in total star particle mass (Msun)", labelpad=2)
plot(t * unit_time_in_cgs / Gyr_in_cgs, (swift_box_star_mass)* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', label='SWIFT')
plot(eagle_time_Gyr[1:], swift_box_star_mass[0] * unit_mass_in_cgs / Msun_in_cgs - eagle_total_mass[:-1],linewidth=0.5,color='k',label='EAGLE test', ls='--')
xlabel("${\\rm Time~[Gyr]}$", labelpad=0)
ylabel("Change in total star particle mass ${[\\rm M_\\odot]}$", labelpad=2)
ticklabel_format(style='sci', axis='y', scilimits=(0,0))
legend()
......@@ -218,36 +220,26 @@ subplot(223)
for j in range(n_elements):
plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (swift_element_mass[1:,j] - swift_element_mass[0,j]) * unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color=colours[j], ms=0.5, label=element_names[j])
plot(eagle_time_Gyr[1:],eagle_total_element_mass[:-1,j],linewidth=1,color=colours[j],linestyle='--')
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in element mass of gas particles (Msun)", labelpad=2)
xlabel("${\\rm Time~[Gyr]}$", labelpad=0)
ylabel("Change in element mass of gas particles ${[\\rm M_\\odot]}$", labelpad=2)
xscale("log")
yscale("log")
legend(bbox_to_anchor=(1.005, 1.), ncol=1, fontsize=8, handlelength=1)
# Box gas metal mass --------------------------------
subplot(224)
plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (swift_box_gas_metal_mass[1:] - swift_box_gas_metal_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', marker = "*", ms=0.5, label='swift')
plot(eagle_time_Gyr[1:],eagle_total_metal_mass[:-1],linewidth=0.5,color='r',label='eagle test')
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in total metal mass of gas particles (Msun)", labelpad=2)
plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (swift_box_gas_metal_mass[1:] - swift_box_gas_metal_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', label='Total')
plot(t * unit_time_in_cgs / Gyr_in_cgs, swift_box_gas_metal_mass_AGB * unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='C0', label='AGB')
plot(t * unit_time_in_cgs / Gyr_in_cgs, swift_box_gas_metal_mass_SNII * unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='C1', label='SNII')
plot(t * unit_time_in_cgs / Gyr_in_cgs, swift_box_gas_metal_mass_SNIa * unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='C2', label='SNIa')
plot(eagle_time_Gyr[1:],eagle_total_metal_mass[:-1],linewidth=0.5,color='k', ls='--')
plot(eagle_time_Gyr[1:],eagle_total_metal_mass_AGB[:-1],linewidth=0.5, color='C0', ls='--')
plot(eagle_time_Gyr[1:],eagle_total_metal_mass_SNII[:-1],linewidth=0.5, color='C1', ls='--')
plot(eagle_time_Gyr[1:],eagle_total_metal_mass_SNIa[:-1],linewidth=0.5, color='C2', ls='--')
legend(loc="center right", ncol=2, fontsize=8)
xlabel("${\\rm Time~[Gyr]}$", labelpad=0)
ylabel("Change in total metal mass of gas particles ${[\\rm M_\\odot]}$", labelpad=2)
ticklabel_format(style='sci', axis='y', scilimits=(0,0))
savefig("box_evolution_Z_%.4f.png"%(Z_star), dpi=200)
# Energy plot
figure()
plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (swift_total_energy[1:] - swift_total_energy[0]) * unit_energy_in_cgs, linewidth=0.5, color='k', label='swift')
plot(eagle_time_Gyr[1:], eagle_energy_SNIa_cgs[:-1], linewidth=0.5, color='b', label='eagle SNIa')
plot(eagle_time_Gyr[1:], eagle_energy_SNII_cgs[:-1], linewidth=0.5, color='c', label='eagle SNII')
plot(eagle_time_Gyr[1:], eagle_energy_ejecta_cgs[:-1], linewidth=0.5, color='y', label='eagle ejecta velocity')
plot(eagle_time_Gyr[1:], eagle_energy_from_mass_cgs[:-1], linewidth=0.5, color='g', label='eagle mass energy')
plot(eagle_time_Gyr[1:], eagle_energy_total_cgs[:-1], linewidth=0.5, color='r', label='eagle total')
plot([0,0], [0,0])
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in internal energy of gas particles (erg)", labelpad=2)
yscale("log")
legend(loc="lower right", ncol=2)
savefig("Energy.png")
examples/SubgridTests/StellarEvolution/run.sh
View file @ 519718d7
......@@ -26,22 +26,22 @@ then
./getSolutions.sh
fi
../../swift --feedback --stars --hydro --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.08 2>&1 | tee output_0p08.log
../../swift --temperature --feedback --stars --hydro --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.08 2>&1 | tee output_0p08.log
python plot_box_evolution.py
../../swift --feedback --stars --hydro --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.04 2>&1 | tee output_0p04.log
../../swift --temperature --feedback --stars --hydro --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.04 2>&1 | tee output_0p04.log
python plot_box_evolution.py
../../swift --feedback --stars --hydro --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.01 2>&1 | tee output_0p01.log
../../swift --temperature --feedback --stars --hydro --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.01 2>&1 | tee output_0p01.log
python plot_box_evolution.py
../../swift --feedback --stars --hydro --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.001 2>&1 | tee output_0p001.log
../../swift --temperature --feedback --stars --hydro --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.001 2>&1 | tee output_0p001.log
python plot_box_evolution.py
../../swift --feedback --stars --hydro --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.0001 2>&1 | tee output_0p0001.log
../../swift --temperature --feedback --stars --hydro --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.0001 2>&1 | tee output_0p0001.log
python plot_box_evolution.py
examples/SubgridTests/StellarEvolution/stellar_evolution.yml
View file @ 519718d7
......@@ -80,7 +80,7 @@ EAGLECooling:
# Properties of the EAGLE feedback and enrichment model.
EAGLEFeedback:
use_SNII_feedback: 0 # Global switch for SNII thermal (stochastic) feedback.
use_SNIa_feedback: 1 # Global switch for SNIa thermal (continuous) feedback.
use_SNIa_feedback: 0 # Global switch for SNIa thermal (continuous) feedback.
use_AGB_enrichment: 1 # Global switch for enrichement from AGB stars.
use_SNII_enrichment: 1 # Global switch for enrichement from SNII stars.
use_SNIa_enrichment: 1 # Global switch for enrichement from SNIa stars.
......
examples/main.c
View file @ 519718d7
......@@ -1050,12 +1050,12 @@ int main(int argc, char *argv[]) {
}
/* Verify that we are not using basic modes incorrectly */
if (with_hydro && N_total[0] == 0) {
if (with_hydro && N_total[swift_type_gas] == 0) {
error(
"ERROR: Running with hydrodynamics but no gas particles found in the "
"ICs!");
}
if (with_gravity && N_total[1] == 0) {
if (with_gravity && N_total[swift_type_count] == 0) {
error(
"ERROR: Running with gravity but no gravity particles found in "
"the ICs!");
......
src/chemistry/EAGLE/chemistry.h
View file @ 519718d7
......@@ -299,6 +299,21 @@ chemistry_get_total_metal_mass_fraction_for_feedback(
return sp->chemistry_data.smoothed_metal_mass_fraction_total;
}
/**
* @brief Returns the abundance array (metal mass fractions) of the
* star particle to be used in feedback/enrichment related routines.
*
* EAGLE uses smooth abundances for everything.
*
* @param sp Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float const*
chemistry_get_metal_mass_fraction_for_feedback(
const struct spart* restrict sp) {
return sp->chemistry_data.smoothed_metal_mass_fraction;
}
/**
* @brief Returns the total metallicity (metal mass fraction) of the
* gas particle to be used in cooling related routines.
......
src/chemistry/none/chemistry.h
View file @ 519718d7
......@@ -174,4 +174,88 @@ __attribute__((always_inline)) INLINE static void chemistry_init_part(
*/
__attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
const struct chemistry_global_data* data, struct spart* restrict sp) {}
/**
* @brief Returns the total metallicity (metal mass fraction) of the
* star particle to be used in feedback/enrichment related routines.
*
* No metallicity treatment here -> return 0.
*
* @param sp Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float
chemistry_get_total_metal_mass_fraction_for_feedback(
const struct spart* restrict sp) {
return 0.f;
}
/**
* @brief Returns the abundance array (metal mass fractions) of the
* star particle to be used in feedback/enrichment related routines.
*
* @param sp Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float const*
chemistry_get_metal_mass_fraction_for_feedback(
const struct spart* restrict sp) {
return NULL;
}
/**
* @brief Returns the total metallicity (metal mass fraction) of the
* gas particle to be used in cooling related routines.
*
* No metallicity treatment here -> return 0.
*
* @param p Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float
chemistry_get_total_metal_mass_fraction_for_cooling(
const struct part* restrict p) {
return 0.f;
}
/**
* @brief Returns the abundance array (metal mass fractions) of the
* gas particle to be used in cooling related routines.
*
* @param p Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float const*
chemistry_get_metal_mass_fraction_for_cooling(const struct part* restrict p) {
return NULL;
}
/**
* @brief Returns the total metallicity (metal mass fraction) of the
* gas particle to be used in star formation related routines.
*
* No metallicity treatment here -> return 0.
*
* @param p Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float
chemistry_get_total_metal_mass_fraction_for_star_formation(
const struct part* restrict p) {
return 0.f;
}
/**
* @brief Returns the abundance array (metal mass fractions) of the
* gas particle to be used in star formation related routines.
*
* @param p Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float const*
chemistry_get_metal_mass_fraction_for_star_formation(
const struct part* restrict p) {
return NULL;
}
#endif /* SWIFT_CHEMISTRY_NONE_H */
src/feedback/EAGLE/feedback.c
View file @ 519718d7
This diff is collapsed.
src/feedback/EAGLE/feedback_iact.h
View file @ 519718d7
......@@ -164,15 +164,11 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
pj->chemistry_data.iron_mass_fraction_from_SNIa =
new_iron_from_SNIa_mass * new_mass_inv;
/* Update mass fraction from SNIa */
const double current_mass_from_SNIa =
pj->chemistry_data.mass_from_SNIa * current_mass;
/* Update mass from SNIa */
const double delta_mass_from_SNIa =
si->feedback_data.to_distribute.mass_from_SNIa * Omega_frac;
const double new_mass_from_SNIa =
current_mass_from_SNIa + delta_mass_from_SNIa;
pj->chemistry_data.mass_from_SNIa = new_mass_from_SNIa * new_mass_inv;
pj->chemistry_data.mass_from_SNIa += delta_mass_from_SNIa;
/* Update metal mass fraction from SNIa */
const double current_metal_mass_from_SNIa =
......@@ -185,15 +181,11 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
pj->chemistry_data.metal_mass_fraction_from_SNIa =
new_metal_mass_from_SNIa * new_mass_inv;
/* Update mass fraction from SNII */
const double current_mass_from_SNII =
pj->chemistry_data.mass_from_SNII * current_mass;
/* Update mass from SNII */
const double delta_mass_from_SNII =
si->feedback_data.to_distribute.mass_from_SNII * Omega_frac;
const double new_mass_from_SNII =
current_mass_from_SNII + delta_mass_from_SNII;
pj->chemistry_data.mass_from_SNII = new_mass_from_SNII * new_mass_inv;
pj->chemistry_data.mass_from_SNII += delta_mass_from_SNII;
/* Update metal mass fraction from SNII */
const double current_metal_mass_from_SNII =
......@@ -206,14 +198,11 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
pj->chemistry_data.metal_mass_fraction_from_SNII =
new_metal_mass_from_SNII * new_mass_inv;
/* Update mass fraction from AGB */
const double current_mass_from_AGB =
pj->chemistry_data.mass_from_AGB * current_mass;
/* Update mass from AGB */
const double delta_mass_from_AGB =
si->feedback_data.to_distribute.mass_from_AGB * Omega_frac;
const double new_mass_from_AGB = current_mass_from_AGB + delta_mass_from_AGB;
pj->chemistry_data.mass_from_AGB = new_mass_from_AGB * new_mass_inv;
pj->chemistry_data.mass_from_AGB += delta_mass_from_AGB;
/* Update metal mass fraction from AGB */
const double current_metal_mass_from_AGB =
......
src/feedback/EAGLE/feedback_properties.h
View file @ 519718d7
......@@ -27,25 +27,28 @@
*/
struct yield_table {
/* Yield table mass bins */
/*! Yield table mass bins */
double *mass;
/* Yield table metallicity bins */
/*! Yield table metallicity (metal mass fractions) bins */
double *metallicity;
/* Array to store yield table resampled by IMF mass bins */
/*! Array to store yield table (individual metals produced by the star)
resampled by IMF mass bins */
double *yield_IMF_resampled;
/* Array to store yield table being read in */
/*! Array to store yield table being read in */
double *yield;
/* Array to store table of ejecta resampled by IMF mass bins */
/*! Array to store table of ejecta (metals alredy in the stars that are
ejected) resampled by IMF mass bins */
double *ejecta_IMF_resampled;
/* Array to store table of ejecta being read in */
/*! Array to store table of ejecta being read in */
double *ejecta;
/* Array to store table of total mass released resampled by IMF mass bins */
/*! Array to store table of total mass released ( metals produced by the star)
resampled by IMF mass bins */
double *total_metals_IMF_resampled;
/* Array to store table of total mass released being read in */
......@@ -133,8 +136,11 @@ struct feedback_props {
/*! Inverse of time-scale of the SNIa decay function in Giga-years */
float SNIa_timescale_Gyr_inv;
/*! Maximal mass used for SNIa feedback (in solar masses) */
double SNIa_max_mass_msun;
/*! Log 10 of the maximal mass used for SNIa feedback (in solar masses) */
float log10_SNIa_max_mass_msun;
double log10_SNIa_max_mass_msun;
/*! Energy released by one supernova type II in cgs units */
double E_SNIa_cgs;
......@@ -165,37 +171,37 @@ struct feedback_props {
/* ------------- Parameters for IMF --------------- */
/*! Array to store calculated IMF */
float *imf;
double *imf;
/*! Arrays to store IMF mass bins */
float *imf_mass_bin;
double *imf_mass_bin;
/*! Arrays to store IMF mass bins (log10)*/
float *imf_mass_bin_log10;
double *imf_mass_bin_log10;
/*! Minimal stellar mass considered by the IMF (in solar masses) */
float imf_min_mass_msun;
double imf_min_mass_msun;
/*! Maximal stellar mass considered by the IMF (in solar masses) */
float imf_max_mass_msun;
double imf_max_mass_msun;
/*! Log 10 of the minimal stellar mass considered by the IMF (in solar masses)
*/
float log10_imf_min_mass_msun;
double log10_imf_min_mass_msun;
/*! Log 10 of the maximal stellar mass considered by the IMF (in solar masses)
*/
float log10_imf_max_mass_msun;
double log10_imf_max_mass_msun;
/* ------------ SNe feedback properties ------------ */
/*! Log 10 of the minimal stellar mass considered for SNII feedback (in solar
* masses) */
float log10_SNII_min_mass_msun;
double log10_SNII_min_mass_msun;
/*! Log 10 of the maximal stellar mass considered for SNII feedback (in solar
* masses) */
float log10_SNII_max_mass_msun;
double log10_SNII_max_mass_msun;
/*! Number of type II supernovae per solar mass */
float num_SNII_per_msun;
......
src/feedback/EAGLE/imf.h
View file @ 519718d7
......@@ -58,7 +58,7 @@ INLINE static void determine_imf_bins(
#endif
const int N_bins = eagle_feedback_N_imf_bins;
const float *imf_bins_log10 = feedback_props->imf_mass_bin_log10;
const double *const imf_bins_log10 = feedback_props->imf_mass_bin_log10;
/* Check whether lower mass is within the IMF mass bin range */
log10_min_mass = max(log10_min_mass, imf_bins_log10[0]);
......@@ -91,19 +91,19 @@ INLINE static void determine_imf_bins(
* yield-weighted integration.
* @param feedback_props the #feedback_props data structure
*/
INLINE static float integrate_imf(const float log10_min_mass,
const float log10_max_mass,
const enum eagle_imf_integration_type mode,
const float *stellar_yields,
const struct feedback_props *feedback_props) {
INLINE static double integrate_imf(
const double log10_min_mass, const double log10_max_mass,
const enum eagle_imf_integration_type mode,
const double stellar_yields[eagle_feedback_N_imf_bins],
const struct feedback_props *feedback_props) {
/* Pull out some common terms */
const float *imf = feedback_props->imf;
const float *imf_mass_bin = feedback_props->imf_mass_bin;
const float *imf_mass_bin_log10 = feedback_props->imf_mass_bin_log10;
const double *imf = feedback_props->imf;
const double *imf_mass_bin = feedback_props->imf_mass_bin;
const double *imf_mass_bin_log10 = feedback_props->imf_mass_bin_log10;
/* IMF mass bin spacing in log10 space. Assumes uniform spacing. */
const float imf_log10_mass_bin_size =
const double imf_log10_mass_bin_size =
imf_mass_bin_log10[1] - imf_mass_bin_log10[0];
/* Determine bins to integrate over based on integration bounds */
......@@ -112,7 +112,7 @@ INLINE static float integrate_imf(const float log10_min_mass,
feedback_props);
/* Array for the integrand */
float integrand[eagle_feedback_N_imf_bins];
double integrand[eagle_feedback_N_imf_bins];
/* Add up the contribution from each of the IMF mass bins */
switch (mode) {
......@@ -153,7 +153,7 @@ INLINE static float integrate_imf(const float log10_min_mass,
}
/* Integrate using trapezoidal rule */
float result = 0.f;
double result = 0.;
for (int i = i_min; i < i_max + 1; i++) {
result += integrand[i];
}
......@@ -163,31 +163,31 @@ INLINE static float integrate_imf(const float log10_min_mass,
result -= 0.5 * (integrand[i_min] + integrand[i_max]);
/* Correct first bin */
const float first_bin_offset =
const double first_bin_offset =
(log10_min_mass - imf_mass_bin_log10[i_min]) / imf_log10_mass_bin_size;