diff --git a/examples/SubgridTests/StellarEvolution/plot_box_evolution.py b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
index aa1f07f98ecede241423bc72a3e0cf7d76963979..eb555bc3c87cfda5f6766ce630254ecafd605d7a 100644
--- a/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
+++ b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
@@ -1,29 +1,30 @@
###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- #
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- #
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- # GNU General Public License for more details.
- #
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program. If not, see <http://www.gnu.org/licenses/>.
- #
- ##############################################################################
-
-# Script used to plot time evolution of gas particle properties. Intended to
+# This file is part of SWIFT.
+# Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+# Script used to plot time evolution of gas particle properties. Intended to
# compare result of feedback due to one star placed in centre of uniform box
# of gas with output from EAGLE feedback test. Can also use as input output
# from SWIFT feedback test (tests/testFeedback) with the appropriate change
# to filepath.
import matplotlib
+
matplotlib.use("Agg")
from pylab import *
from scipy import stats
@@ -33,31 +34,34 @@ import glob
import os.path
# Plot parameters
-params = {'axes.labelsize': 10,
-'axes.titlesize': 10,
-'font.size': 12,
-'legend.fontsize': 12,
-'xtick.labelsize': 10,
-'ytick.labelsize': 10,
-'text.usetex': True,
- 'figure.figsize' : (9.90,6.45),
-'figure.subplot.left' : 0.05,
-'figure.subplot.right' : 0.995,
-'figure.subplot.bottom' : 0.06,
-'figure.subplot.top' : 0.92,
-'figure.subplot.wspace' : 0.25,
-'figure.subplot.hspace' : 0.2,
-'lines.markersize' : 6,
-'lines.linewidth' : 3.,
-'text.latex.unicode': True
+params = {
+ "axes.labelsize": 10,
+ "axes.titlesize": 10,
+ "font.size": 12,
+ "legend.fontsize": 12,
+ "xtick.labelsize": 10,
+ "ytick.labelsize": 10,
+ "text.usetex": True,
+ "figure.figsize": (9.90, 6.45),
+ "figure.subplot.left": 0.05,
+ "figure.subplot.right": 0.995,
+ "figure.subplot.bottom": 0.06,
+ "figure.subplot.top": 0.92,
+ "figure.subplot.wspace": 0.25,
+ "figure.subplot.hspace": 0.2,
+ "lines.markersize": 6,
+ "lines.linewidth": 3.0,
+ "text.latex.unicode": True,
}
rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
# Number of snapshots and elements
-newest_snap_name = max(glob.glob('stellar_evolution_*.hdf5'))#, key=os.path.getctime)
-n_snapshots = int(newest_snap_name.replace('stellar_evolution_','').replace('.hdf5','')) + 1
+newest_snap_name = max(glob.glob("stellar_evolution_*.hdf5")) # , key=os.path.getctime)
+n_snapshots = (
+ int(newest_snap_name.replace("stellar_evolution_", "").replace(".hdf5", "")) + 1
+)
n_elements = 9
# Read the simulation data
@@ -72,7 +76,9 @@ git = sim["Code"].attrs["Git Revision"]
stellar_mass = sim["/PartType4/Masses"][0]
E_SNII_cgs = double(sim["/Parameters"].attrs["EAGLEFeedback:SNII_energy_erg"])
E_SNIa_cgs = double(sim["/Parameters"].attrs["EAGLEFeedback:SNIa_energy_erg"])
-ejecta_vel_cgs = double(sim["/Parameters"].attrs["EAGLEFeedback:AGB_ejecta_velocity_km_p_s"]) * 1e5
+ejecta_vel_cgs = (
+ double(sim["/Parameters"].attrs["EAGLEFeedback:AGB_ejecta_velocity_km_p_s"]) * 1e5
+)
# Units
unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
@@ -81,11 +87,11 @@ unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
unit_temp_in_cgs = sim["/Units"].attrs["Unit temperature in cgs (U_T)"]
unit_vel_in_cgs = unit_length_in_cgs / unit_time_in_cgs
unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
-unit_density_in_cgs = unit_mass_in_cgs*unit_length_in_cgs**-3
-unit_pressure_in_cgs = unit_mass_in_cgs/unit_length_in_cgs*unit_time_in_cgs**-2
-unit_int_energy_in_cgs = unit_energy_in_cgs/unit_mass_in_cgs
-unit_entropy_in_cgs = unit_energy_in_cgs/unit_temp_in_cgs
-Gyr_in_cgs = 1e9 * 365. * 24 * 3600.
+unit_density_in_cgs = unit_mass_in_cgs * unit_length_in_cgs ** -3
+unit_pressure_in_cgs = unit_mass_in_cgs / unit_length_in_cgs * unit_time_in_cgs ** -2
+unit_int_energy_in_cgs = unit_energy_in_cgs / unit_mass_in_cgs
+unit_entropy_in_cgs = unit_energy_in_cgs / unit_temp_in_cgs
+Gyr_in_cgs = 1e9 * 365.0 * 24 * 3600.0
Msun_in_cgs = 1.98848e33
# Declare arrays to store SWIFT data
@@ -98,148 +104,247 @@ 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_element_mass = zeros((n_snapshots, n_elements))
swift_internal_energy = zeros(n_snapshots)
swift_kinetic_energy = zeros(n_snapshots)
swift_total_energy = zeros(n_snapshots)
-swift_mean_u_start = 0.
+swift_mean_u_start = 0.0
t = zeros(n_snapshots)
# Read data from snapshots
for i in range(n_snapshots):
- #print("reading snapshot "+str(i))
- sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
- t[i] = sim["/Header"].attrs["Time"][0]
+ # print("reading snapshot "+str(i))
+ sim = h5py.File("stellar_evolution_%04d.hdf5" % i, "r")
+ t[i] = sim["/Header"].attrs["Time"][0]
+
+ masses = sim["/PartType0/Masses"][:]
+ swift_box_gas_mass[i] = np.sum(masses)
- 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"][:]
- 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)
- 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)
- Z_star = sim["/PartType4/MetalMassFractions"][0]
- star_masses = sim["/PartType4/Masses"][:]
- swift_box_star_mass[i] = np.sum(star_masses)
+ metallicities = sim["/PartType0/MetalMassFractions"][:]
+ swift_box_gas_metal_mass[i] = np.sum(metallicities * masses)
- 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"][:]
- 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)
- 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)
- element_abundances = sim["/PartType0/ElementMassFractions"][:][:]
- for j in range(n_elements):
- swift_element_mass[i,j] = np.sum(element_abundances[:,j] * masses)
+ v = sim["/PartType0/Velocities"][:, :]
+ v2 = v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2
+ u = sim["/PartType0/InternalEnergies"][:]
+ swift_internal_energy[i] = np.sum(masses * u)
+ swift_kinetic_energy[i] = np.sum(0.5 * masses * v2)
+ swift_total_energy[i] = swift_kinetic_energy[i] + swift_internal_energy[i]
- v = sim["/PartType0/Velocities"][:,:]
- v2 = v[:,0]**2 + v[:,1]**2 + v[:,2]**2
- u = sim["/PartType0/InternalEnergies"][:]
- swift_internal_energy[i] = np.sum(masses * u)
- swift_kinetic_energy[i] = np.sum(0.5 * masses * v2)
- swift_total_energy[i] = swift_kinetic_energy[i] + swift_internal_energy[i]
+ if i == 0:
+ swift_mean_u_start = np.mean(u)
- if i == 0:
- swift_mean_u_start = np.mean(u)
-
- sim.close()
+ sim.close()
# Read expected yields from EAGLE. Choose which file to use based on metallicity used when
# running SWIFT (can be specified in yml file)
-filename = "./StellarEvolutionSolution/Z_%.4f/StellarEvolutionTotal.txt"%Z_star
+filename = "./StellarEvolutionSolution/Z_%.4f/StellarEvolutionTotal.txt" % Z_star
# Read EAGLE test output
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
+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
+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
+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
+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
+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
-filename = "./StellarEvolutionSolution/Z_%.4f/StellarEvolutionIa.txt"%Z_star
+filename = "./StellarEvolutionSolution/Z_%.4f/StellarEvolutionIa.txt" % Z_star
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
+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()
-suptitle("Star metallicity Z = %.4f"%Z_star)
+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', 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='--')
+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))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
# 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', 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='--')
+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))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
legend()
# Box gas element mass --------------------------------
-colours = ['k','r','g','b','c','y','m','skyblue','plum']
-element_names = ['H','He','C','N','O','Ne','Mg','Si','Fe']
+colours = ["k", "r", "g", "b", "c", "y", "m", "skyblue", "plum"]
+element_names = ["H", "He", "C", "N", "O", "Ne", "Mg", "Si", "Fe"]
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='--')
+ 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 ${[\\rm M_\\odot]}$", labelpad=2)
xscale("log")
yscale("log")
-legend(bbox_to_anchor=(1.005, 1.), ncol=1, fontsize=8, handlelength=1)
+legend(bbox_to_anchor=(1.005, 1.0), 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', 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='--')
+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)
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
+savefig("box_evolution_Z_%.4f.png" % (Z_star), dpi=200)
diff --git a/examples/SubgridTests/StellarEvolution/run.sh b/examples/SubgridTests/StellarEvolution/run.sh
index 032083db5055dadac8c1bdcdc6cad3e2e03498a1..28509b8fe52a39e2a2d4049ee62813747a72b0b1 100755
--- a/examples/SubgridTests/StellarEvolution/run.sh
+++ b/examples/SubgridTests/StellarEvolution/run.sh
@@ -28,20 +28,20 @@ fi
../../swift --temperature --feedback --stars --hydro --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.08 -P EAGLEChemistry:init_abundance_Hydrogen:0.71 -P EAGLEChemistry:init_abundance_Helium:0.21 2>&1 | tee output_0p08.log
-python plot_box_evolution.py
+python3 plot_box_evolution.py
../../swift --temperature --feedback --stars --hydro --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.04 -P EAGLEChemistry:init_abundance_Hydrogen:0.74 -P EAGLEChemistry:init_abundance_Helium:0.23 2>&1 | tee output_0p04.log
-python plot_box_evolution.py
+python3 plot_box_evolution.py
../../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
+python3 plot_box_evolution.py
../../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
+python3 plot_box_evolution.py
../../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
+python3 plot_box_evolution.py
diff --git a/examples/SubgridTests/StellarEvolution/stellar_evolution.yml b/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
index 9836aa5931dddb5d7f8e5036d246ab9ec51e70b1..86bc8e234ebba2af6581f74c90538a539c48126e 100644
--- a/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
+++ b/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
@@ -90,6 +90,7 @@ EAGLEFeedback:
IMF_max_mass_Msun: 100.0 # Maximal stellar mass considered for the Chabrier IMF in solar masses.
SNII_min_mass_Msun: 6.0 # Minimal mass considered for SNII stars in solar masses.
SNII_max_mass_Msun: 100.0 # Maximal mass considered for SNII stars in solar masses.
+ SNII_feedback_model: Random # Feedback modes: Random, Isotropic, MinimumDistance, MinimumDensity
SNII_sampled_delay: 0 # Sample the SNII lifetimes to do feedback.
SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event when not sampling.
SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
@@ -100,7 +101,8 @@ EAGLEFeedback:
SNII_energy_fraction_n_0_H_p_cm3: 0.67 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
SNII_energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction.
SNII_energy_fraction_n_n: 0.8686 # Power-law for the birth density dependance of the SNII energy fraction.
- SNII_energy_fraction_use_birth_props: 0 # Are we using the density and metallicity at birth to compute f_E or at feedback time?
+ SNII_energy_fraction_use_birth_density: 0 # Are we using the density and metallicity at birth to compute f_E or at feedback time?
+ SNII_energy_fraction_use_birth_metallicity: 0 # Are we using the density and metallicity at birth to compute f_E or at feedback time?
SNIa_DTD: Exponential # Use the EAGLE-Ref SNIa DTD.
SNIa_DTD_delay_Gyr: 0.04 # Age of the most massive SNIa in Gyr.
SNIa_DTD_exp_timescale_Gyr: 2.0 # Time-scale of the exponential decay of the SNIa rates in Gyr.
diff --git a/src/cooling/grackle/cooling.c b/src/cooling/grackle/cooling.c
index 85b931f5e049b8cac220cf68f20c32f4a18186d4..84b9fa4457f5cce15ea0cdb61f530d8548c43ae0 100644
--- a/src/cooling/grackle/cooling.c
+++ b/src/cooling/grackle/cooling.c
@@ -256,34 +256,6 @@ void cooling_first_init_part(const struct phys_const* restrict phys_const,
#endif
}
-/**
- * @brief Returns the subgrid temperature of a particle.
- *
- * This model has no subgrid quantity. We return an error.
- *
- * @param p The particle.
- * @param xp The extended particle data.
- */
-float cooling_get_subgrid_temperature(const struct part* p,
- const struct xpart* xp) {
- error("This cooling model does not use subgrid quantities!");
- return -1.f;
-}
-
-/**
- * @brief Returns the subgrid density of a particle.
- *
- * This model has no subgrid quantity. We return an error.
- *
- * @param p The particle.
- * @param xp The extended particle data.
- */
-float cooling_get_subgrid_density(const struct part* p,
- const struct xpart* xp) {
- error("This cooling model does not use subgrid quantities!");
- return -1.f;
-}
-
/**
* @brief Returns the total radiated energy by this particle.
*
diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h
index 4c58ef8fb1b2cddcfd965cf7a6c90aed2cbbe82d..22dd1d642c6d581c821cbbb8ac5bb33feacce7fa 100644
--- a/src/cooling/grackle/cooling.h
+++ b/src/cooling/grackle/cooling.h
@@ -68,10 +68,33 @@ void cooling_first_init_part(const struct phys_const* restrict phys_const,
const struct part* restrict p,
struct xpart* restrict xp);
-float cooling_get_subgrid_temperature(const struct part* p,
- const struct xpart* xp);
+/**
+ * @brief Returns the subgrid temperature of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+INLINE static float cooling_get_subgrid_temperature(const struct part* p,
+ const struct xpart* xp) {
+ error("This cooling model does not use subgrid quantities!");
+ return -1.f;
+}
-float cooling_get_subgrid_density(const struct part* p, const struct xpart* xp);
+/**
+ * @brief Returns the subgrid density of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+INLINE static float cooling_get_subgrid_density(const struct part* p,
+ const struct xpart* xp) {
+ error("This cooling model does not use subgrid quantities!");
+ return -1.f;
+}
float cooling_get_radiated_energy(const struct xpart* restrict xp);
void cooling_print_backend(const struct cooling_function_data* cooling);
diff --git a/src/feedback/EAGLE/feedback.h b/src/feedback/EAGLE/feedback.h
index 51294394e625f88a0cf41cafbbf83a1634b537fd..b8a905cac1c8100ba2c748c000903df7059d84ff 100644
--- a/src/feedback/EAGLE/feedback.h
+++ b/src/feedback/EAGLE/feedback.h
@@ -170,10 +170,9 @@ __attribute__((always_inline)) INLINE static void feedback_prepare_spart(
struct spart* sp, const struct feedback_props* feedback_props) {}
/**
- * @brief Evolve the stellar properties of a #spart.
+ * @brief Prepare a #spart for the feedback task.
*
- * This function allows for example to compute the SN rate before sending
- * this information to a different MPI rank.
+ * In EAGLE, this function evolves the stellar properties of a #spart.
*
* @param sp The particle to act upon
* @param feedback_props The #feedback_props structure.
@@ -187,7 +186,7 @@ __attribute__((always_inline)) INLINE static void feedback_prepare_spart(
* @param ti_begin The integer time at the beginning of the step.
* @param with_cosmology Are we running with cosmology on?
*/
-__attribute__((always_inline)) INLINE static void feedback_evolve_spart(
+__attribute__((always_inline)) INLINE static void feedback_prepare_feedback(
struct spart* restrict sp, const struct feedback_props* feedback_props,
const struct cosmology* cosmo, const struct unit_system* us,
const struct phys_const* phys_const, const double star_age_beg_step,
@@ -225,16 +224,27 @@ __attribute__((always_inline)) INLINE static void feedback_evolve_spart(
/**
* @brief Will this star particle want to do feedback during the next time-step?
*
- * @param sp The star of interest.
- * @param feedback_props The properties of the feedback model.
- * @param with_cosmology Are we running a cosmological problem?
- * @param cosmo The cosmological model.
- * @param time The current time (since the start of the run / Big Bang).
+ * This is called in the time step task and increases counters of time-steps
+ * that have been performed.
+ *
+ * @param sp The particle to act upon
+ * @param feedback_props The #feedback_props structure.
+ * @param cosmo The current cosmological model.
+ * @param us The unit system.
+ * @param phys_const The #phys_const.
+ * @param star_age_beg_step The age of the star at the star of the time-step in
+ * internal units.
+ * @param dt The time-step size of this star in internal units.
+ * @param time The physical time in internal units.
+ * @param ti_begin The integer time at the beginning of the step.
+ * @param with_cosmology Are we running with cosmology on?
*/
-__attribute__((always_inline)) INLINE static int feedback_will_do_feedback(
- struct spart* restrict sp, const struct feedback_props* feedback_props,
- const int with_cosmology, const struct cosmology* cosmo,
- const double time) {
+__attribute__((always_inline)) INLINE static void feedback_will_do_feedback(
+ struct spart* sp, const struct feedback_props* feedback_props,
+ const int with_cosmology, const struct cosmology* cosmo, const double time,
+ const struct unit_system* us, const struct phys_const* phys_const,
+ const double star_age_beg_step, const double dt,
+ const integertime_t ti_begin) {
/* Special case for new-born stars */
if (with_cosmology) {
@@ -242,18 +252,12 @@ __attribute__((always_inline)) INLINE static int feedback_will_do_feedback(
/* Set the counter to "let's do enrichment" */
sp->count_since_last_enrichment = 0;
-
- /* Say we want to do feedback */
- return 1;
}
} else {
if (sp->birth_time == (float)time) {
/* Set the counter to "let's do enrichment" */
sp->count_since_last_enrichment = 0;
-
- /* Say we want to do feedback */
- return 1;
}
}
@@ -272,9 +276,6 @@ __attribute__((always_inline)) INLINE static int feedback_will_do_feedback(
/* Set the counter to "let's do enrichment" */
sp->count_since_last_enrichment = 0;
- /* Say we want to do feedback */
- return 1;
-
} else {
/* Increment counter */
@@ -285,14 +286,6 @@ __attribute__((always_inline)) INLINE static int feedback_will_do_feedback(
/* Reset counter */
sp->count_since_last_enrichment = 0;
-
- /* Say we want to do feedback */
- return 1;
-
- } else {
-
- /* Say we don't want to do feedback */
- return 0;
}
}
}
diff --git a/src/feedback/GEAR/feedback.c b/src/feedback/GEAR/feedback.c
index c418465ab5097866094adc85216c5e500ad00b21..38662719854f7a507173ed4bef4e0f9620939e48 100644
--- a/src/feedback/GEAR/feedback.c
+++ b/src/feedback/GEAR/feedback.c
@@ -87,21 +87,69 @@ void feedback_update_part(struct part* restrict p, struct xpart* restrict xp,
}
/**
- * @brief Should we do feedback for this star?
+ * @brief Will this star particle want to do feedback during the next time-step?
*
- * @param sp The star to consider.
- * @param feedback_props The #feedback_props.
- * @param with_cosmology Is the cosmology switch on?
- * @param cosmo The #cosmology.
- * @param time The current time.
+ * This is called in the time step task.
+ *
+ * In GEAR, we compute the full stellar evolution here.
+ *
+ * @param sp The particle to act upon
+ * @param feedback_props The #feedback_props structure.
+ * @param cosmo The current cosmological model.
+ * @param us The unit system.
+ * @param phys_const The #phys_const.
+ * @param star_age_beg_step The age of the star at the star of the time-step in
+ * internal units.
+ * @param dt The time-step size of this star in internal units.
+ * @param time The physical time in internal units.
+ * @param ti_begin The integer time at the beginning of the step.
+ * @param with_cosmology Are we running with cosmology on?
*/
-int feedback_will_do_feedback(const struct spart* sp,
- const struct feedback_props* feedback_props,
- const int with_cosmology,
- const struct cosmology* cosmo,
- const double time) {
+void feedback_will_do_feedback(struct spart* sp,
+ const struct feedback_props* feedback_props,
+ const int with_cosmology,
+ const struct cosmology* cosmo, const double time,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ const double star_age_beg_step, const double dt,
+ const integertime_t ti_begin) {
+
+ /* Zero the energy of supernovae */
+ sp->feedback_data.energy_ejected = 0;
+ sp->feedback_data.will_do_feedback = 0;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (sp->birth_time == -1.) error("Evolving a star particle that should not!");
+
+ if (star_age_beg_step < -1e-6) {
+ error("Negative age for a star");
+ }
+#endif
+ /* Has this star been around for a while ? */
+ if (star_age_beg_step + dt <= 0.) return;
- return (sp->birth_time != -1.);
+ const double star_age_beg_step_safe =
+ star_age_beg_step < 0 ? 0 : star_age_beg_step;
+
+ /* Pick the correct table. (if only one table, threshold is < 0) */
+ const float metal =
+ chemistry_get_star_total_metal_mass_fraction_for_feedback(sp);
+ const float threshold = feedback_props->metallicity_max_first_stars;
+
+ const struct stellar_model* model =
+ metal < threshold ? &feedback_props->stellar_model_first_stars
+ : &feedback_props->stellar_model;
+
+ /* Compute the stellar evolution */
+ stellar_evolution_evolve_spart(sp, model, cosmo, us, phys_const, ti_begin,
+ star_age_beg_step_safe, dt);
+
+ /* Transform the number of SN to the energy */
+ sp->feedback_data.energy_ejected =
+ sp->feedback_data.number_sn * feedback_props->energy_per_supernovae;
+
+ /* Set the particle as doing some feedback */
+ sp->feedback_data.will_do_feedback = sp->feedback_data.energy_ejected != 0.;
}
/**
@@ -116,14 +164,7 @@ int feedback_is_active(const struct spart* sp, const double time,
const struct cosmology* cosmo,
const int with_cosmology) {
- // TODO improve this with estimates for SNII and SNIa
- if (sp->birth_time == -1.) return 0;
-
- if (with_cosmology) {
- return ((double)cosmo->a) > sp->birth_scale_factor;
- } else {
- return time > sp->birth_time;
- }
+ return sp->feedback_data.will_do_feedback;
}
/**
@@ -157,16 +198,31 @@ void feedback_init_spart(struct spart* sp) {
}
/**
- * @brief Prepares a star's feedback field before computing what
- * needs to be distributed.
+ * @brief Reset the feedback field when the spart is not
+ * in a correct state for feeedback_will_do_feedback.
+ *
+ * This function is called in the timestep task.
*/
-void feedback_reset_feedback(struct spart* sp,
- const struct feedback_props* feedback_props) {
+void feedback_init_after_star_formation(
+ struct spart* sp, const struct feedback_props* feedback_props) {
+ feedback_init_spart(sp);
/* Zero the energy of supernovae */
sp->feedback_data.energy_ejected = 0;
+
+ /* Activate the feedback loop for the first step */
+ sp->feedback_data.will_do_feedback = 1;
}
+/**
+ * @brief Prepares a star's feedback field before computing what
+ * needs to be distributed.
+ *
+ * This is called in the stars ghost.
+ */
+void feedback_reset_feedback(struct spart* sp,
+ const struct feedback_props* feedback_props) {}
+
/**
* @brief Initialises the s-particles feedback props for the first time
*
@@ -181,7 +237,11 @@ void feedback_first_init_spart(struct spart* sp,
feedback_init_spart(sp);
- feedback_reset_feedback(sp, feedback_props);
+ /* Zero the energy of supernovae */
+ sp->feedback_data.energy_ejected = 0;
+
+ /* Activate the feedback loop for the first step */
+ sp->feedback_data.will_do_feedback = 1;
}
/**
@@ -197,10 +257,11 @@ void feedback_prepare_spart(struct spart* sp,
const struct feedback_props* feedback_props) {}
/**
- * @brief Evolve the stellar properties of a #spart.
+ * @brief Prepare a #spart for the feedback task.
*
- * This function compute the SN rate and yields before sending
- * this information to a different MPI rank.
+ * This is called in the stars ghost task.
+ *
+ * In here, we only need to add the missing coefficients.
*
* @param sp The particle to act upon
* @param feedback_props The #feedback_props structure.
@@ -214,50 +275,18 @@ void feedback_prepare_spart(struct spart* sp,
* @param ti_begin The integer time at the beginning of the step.
* @param with_cosmology Are we running with cosmology on?
*/
-void feedback_evolve_spart(struct spart* restrict sp,
- const struct feedback_props* feedback_props,
- const struct cosmology* cosmo,
- const struct unit_system* us,
- const struct phys_const* phys_const,
- const double star_age_beg_step, const double dt,
- const double time, const integertime_t ti_begin,
- const int with_cosmology) {
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (sp->birth_time == -1.) error("Evolving a star particle that should not!");
-
- if (star_age_beg_step < -1e-6) {
- error("Negative age for a star");
- }
-#endif
- const double star_age_beg_step_safe =
- star_age_beg_step < 0 ? 0 : star_age_beg_step;
-
- /* Reset the feedback */
- feedback_reset_feedback(sp, feedback_props);
-
+void feedback_prepare_feedback(struct spart* restrict sp,
+ const struct feedback_props* feedback_props,
+ const struct cosmology* cosmo,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ const double star_age_beg_step, const double dt,
+ const double time, const integertime_t ti_begin,
+ const int with_cosmology) {
/* Add missing h factor */
const float hi_inv = 1.f / sp->h;
const float hi_inv_dim = pow_dimension(hi_inv); /* 1/h^d */
-
sp->feedback_data.enrichment_weight *= hi_inv_dim;
-
- /* Pick the correct table. (if only one table, threshold is < 0) */
- const float metal =
- chemistry_get_star_total_metal_mass_fraction_for_feedback(sp);
- const float threshold = feedback_props->metallicity_max_first_stars;
-
- const struct stellar_model* model =
- metal < threshold ? &feedback_props->stellar_model_first_stars
- : &feedback_props->stellar_model;
-
- /* Compute the stellar evolution */
- stellar_evolution_evolve_spart(sp, model, cosmo, us, phys_const, ti_begin,
- star_age_beg_step_safe, dt);
-
- /* Transform the number of SN to the energy */
- sp->feedback_data.energy_ejected =
- sp->feedback_data.number_sn * feedback_props->energy_per_supernovae;
}
/**
diff --git a/src/feedback/GEAR/feedback.h b/src/feedback/GEAR/feedback.h
index c4b26b29e746fc192ec82e02fb077faa09ab8611..bd0306b158e41ddde1b4d825cc5a30fa81154205 100644
--- a/src/feedback/GEAR/feedback.h
+++ b/src/feedback/GEAR/feedback.h
@@ -32,10 +32,14 @@
void feedback_update_part(struct part* restrict p, struct xpart* restrict xp,
const struct engine* restrict e);
-int feedback_will_do_feedback(const struct spart* sp,
- const struct feedback_props* feedback_props,
- const int with_cosmology,
- const struct cosmology* cosmo, const double time);
+void feedback_will_do_feedback(struct spart* sp,
+ const struct feedback_props* feedback_props,
+ const int with_cosmology,
+ const struct cosmology* cosmo, const double time,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ const double star_age_beg_step, const double dt,
+ const integertime_t ti_begin);
int feedback_is_active(const struct spart* sp, const double time,
const struct cosmology* cosmo, const int with_cosmology);
@@ -46,20 +50,22 @@ double feedback_get_enrichment_timestep(const struct spart* sp,
const double dt_star);
void feedback_init_spart(struct spart* sp);
+void feedback_init_after_star_formation(
+ struct spart* sp, const struct feedback_props* feedback_props);
void feedback_reset_feedback(struct spart* sp,
const struct feedback_props* feedback_props);
void feedback_first_init_spart(struct spart* sp,
const struct feedback_props* feedback_props);
void feedback_prepare_spart(struct spart* sp,
const struct feedback_props* feedback_props);
-void feedback_evolve_spart(struct spart* restrict sp,
- const struct feedback_props* feedback_props,
- const struct cosmology* cosmo,
- const struct unit_system* us,
- const struct phys_const* phys_const,
- const double star_age_beg_step, const double dt,
- const double time, const integertime_t ti_begin,
- const int with_cosmology);
+void feedback_prepare_feedback(struct spart* restrict sp,
+ const struct feedback_props* feedback_props,
+ const struct cosmology* cosmo,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ const double star_age_beg_step, const double dt,
+ const double time, const integertime_t ti_begin,
+ const int with_cosmology);
void feedback_struct_dump(const struct feedback_props* feedback, FILE* stream);
void feedback_struct_restore(struct feedback_props* feedback, FILE* stream);
void feedback_clean(struct feedback_props* feedback);
diff --git a/src/feedback/GEAR/feedback_struct.h b/src/feedback/GEAR/feedback_struct.h
index bc000bea0dc0f45be9cec35b928742f9f4ae64f1..f047462555c04a4ee1a0e4bd3a45284e22cf61fa 100644
--- a/src/feedback/GEAR/feedback_struct.h
+++ b/src/feedback/GEAR/feedback_struct.h
@@ -58,6 +58,9 @@ struct feedback_spart_data {
/*! Chemical composition of the mass ejected */
double metal_mass_ejected[GEAR_CHEMISTRY_ELEMENT_COUNT];
+
+ /*! Does the particle needs the feedback loop? */
+ char will_do_feedback;
};
#endif /* SWIFT_FEEDBACK_STRUCT_GEAR_H */
diff --git a/src/feedback/none/feedback.h b/src/feedback/none/feedback.h
index e5d42f456b4b1a9246f5e24c22aeb3136c9fc799..1a0e869417fe5fc2c1d691af7dd9d811310bf95b 100644
--- a/src/feedback/none/feedback.h
+++ b/src/feedback/none/feedback.h
@@ -104,6 +104,8 @@ INLINE static double feedback_get_enrichment_timestep(
/**
* @brief Prepares a star's feedback field before computing what
* needs to be distributed.
+ *
+ * This is called in the stars ghost.
*/
__attribute__((always_inline)) INLINE static void feedback_reset_feedback(
struct spart* sp, const struct feedback_props* feedback_props) {}
@@ -133,8 +135,9 @@ __attribute__((always_inline)) INLINE static void feedback_prepare_spart(
struct spart* sp, const struct feedback_props* feedback_props) {}
/**
- * @brief Evolve the stellar properties of a #spart.
+ * @brief Prepare a #spart for the feedback task.
*
+ * This is called in the stars ghost task.
* This function allows for example to compute the SN rate before sending
* this information to a different MPI rank.
*
@@ -150,7 +153,7 @@ __attribute__((always_inline)) INLINE static void feedback_prepare_spart(
* @param ti_begin The integer time at the beginning of the step.
* @param with_cosmology Are we running with cosmology on?
*/
-__attribute__((always_inline)) INLINE static void feedback_evolve_spart(
+__attribute__((always_inline)) INLINE static void feedback_prepare_feedback(
struct spart* restrict sp, const struct feedback_props* feedback_props,
const struct cosmology* cosmo, const struct unit_system* us,
const struct phys_const* phys_const, const double star_age_beg_step,
@@ -160,18 +163,26 @@ __attribute__((always_inline)) INLINE static void feedback_evolve_spart(
/**
* @brief Will this star particle want to do feedback during the next time-step?
*
- * @param sp The star of interest.
- * @param feedback_props The properties of the feedback model.
- * @param with_cosmology Are we running with cosmological time integration?
- * @param cosmo The #cosmology object.
- * @param time The current time (since the Big Bang).
+ * This is called in the time step task.
+ *
+ * @param sp The particle to act upon
+ * @param feedback_props The #feedback_props structure.
+ * @param cosmo The current cosmological model.
+ * @param us The unit system.
+ * @param phys_const The #phys_const.
+ * @param star_age_beg_step The age of the star at the star of the time-step in
+ * internal units.
+ * @param dt The time-step size of this star in internal units.
+ * @param time The physical time in internal units.
+ * @param ti_begin The integer time at the beginning of the step.
+ * @param with_cosmology Are we running with cosmology on?
*/
-__attribute__((always_inline)) INLINE static int feedback_will_do_feedback(
- struct spart* restrict sp, const struct feedback_props* feedback_props,
- const int with_cosmology, const struct cosmology* cosmo,
- const double time) {
- return 1;
-}
+__attribute__((always_inline)) INLINE static void feedback_will_do_feedback(
+ const struct spart* sp, const struct feedback_props* feedback_props,
+ const int with_cosmology, const struct cosmology* cosmo, const double time,
+ const struct unit_system* us, const struct phys_const* phys_const,
+ const double star_age_beg_step, const double dt,
+ const integertime_t ti_begin) {}
/**
* @brief Clean-up the memory allocated for the feedback routines
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index 8ef7214d80f375cfdd6ace1fcd521c085231e256..b98e47a0f31c4c5f235f6f20589597015934fb7e 100644
--- a/src/runner_ghost.c
+++ b/src/runner_ghost.c
@@ -242,9 +242,10 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
star_age_end_of_step - dt_enrichment;
/* Compute the stellar evolution */
- feedback_evolve_spart(sp, feedback_props, cosmo, us, phys_const,
- star_age_beg_of_step, dt_enrichment,
- e->time, ti_begin, with_cosmology);
+ feedback_prepare_feedback(sp, feedback_props, cosmo, us,
+ phys_const, star_age_beg_of_step,
+ dt_enrichment, e->time, ti_begin,
+ with_cosmology);
} else {
/* Reset the feedback fields of the star particle */
@@ -386,9 +387,9 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
star_age_end_of_step - dt_enrichment;
/* Compute the stellar evolution */
- feedback_evolve_spart(sp, feedback_props, cosmo, us, phys_const,
- star_age_beg_of_step, dt_enrichment, e->time,
- ti_begin, with_cosmology);
+ feedback_prepare_feedback(sp, feedback_props, cosmo, us, phys_const,
+ star_age_beg_of_step, dt_enrichment,
+ e->time, ti_begin, with_cosmology);
} else {
/* Reset the feedback fields of the star particle */
diff --git a/src/runner_time_integration.c b/src/runner_time_integration.c
index db7f4a23ba743683137308680381ff86cc5608ba..b7ab9c68793bf79c6e0317d265faf8d5e87f1ce3 100644
--- a/src/runner_time_integration.c
+++ b/src/runner_time_integration.c
@@ -624,6 +624,10 @@ void runner_do_timestep(struct runner *r, struct cell *c, const int timer) {
const struct engine *e = r->e;
const integertime_t ti_current = e->ti_current;
+ const struct cosmology *cosmo = e->cosmology;
+ const struct feedback_props *feedback_props = e->feedback_props;
+ const struct unit_system *us = e->internal_units;
+ const struct phys_const *phys_const = e->physical_constants;
const int with_cosmology = (e->policy & engine_policy_cosmology);
const int with_feedback = (e->policy & engine_policy_feedback);
const int count = c->hydro.count;
@@ -839,9 +843,40 @@ void runner_do_timestep(struct runner *r, struct cell *c, const int timer) {
sp->gpart->time_bin = get_time_bin(ti_new_step);
/* Update feedback related counters */
- if (with_feedback)
- feedback_will_do_feedback(sp, e->feedback_props, with_cosmology,
- e->cosmology, e->time);
+ if (with_feedback) {
+ const integertime_t ti_step = get_integer_timestep(sp->time_bin);
+ const integertime_t ti_begin_star =
+ get_integer_time_begin(e->ti_current, sp->time_bin);
+
+ /* Get particle time-step */
+ double dt_star;
+ if (with_cosmology) {
+ dt_star = cosmology_get_delta_time(e->cosmology, ti_begin_star,
+ ti_begin_star + ti_step);
+ } else {
+ dt_star = get_timestep(sp->time_bin, e->time_base);
+ }
+
+ /* Calculate age of the star at current time */
+ double star_age_end_of_step;
+ if (with_cosmology) {
+ star_age_end_of_step = cosmology_get_delta_time_from_scale_factors(
+ cosmo, (double)sp->birth_scale_factor, cosmo->a);
+ } else {
+ star_age_end_of_step = e->time - (double)sp->birth_time;
+ }
+
+ /* Get the length of the enrichment time-step */
+ const double dt_enrichment = feedback_get_enrichment_timestep(
+ sp, with_cosmology, cosmo, e->time, dt_star);
+ const double star_age_beg_of_step =
+ star_age_end_of_step - dt_enrichment;
+
+ /* Compute the stellar evolution */
+ feedback_will_do_feedback(
+ sp, feedback_props, with_cosmology, cosmo, e->time, us,
+ phys_const, star_age_beg_of_step, dt_enrichment, ti_begin_star);
+ }
/* Number of updated s-particles */
s_updated++;
diff --git a/src/star_formation/GEAR/star_formation.h b/src/star_formation/GEAR/star_formation.h
index 7db0748ca7eaa3ae119eeeecff9c26208e694d80..c69c7a962dd246289cb8ba26df5858d10d4497b3 100644
--- a/src/star_formation/GEAR/star_formation.h
+++ b/src/star_formation/GEAR/star_formation.h
@@ -27,6 +27,7 @@
#include "engine.h"
#include "entropy_floor.h"
#include "error.h"
+#include "feedback.h"
#include "hydro_properties.h"
#include "parser.h"
#include "part.h"
@@ -307,6 +308,9 @@ INLINE static void star_formation_copy_properties(
const struct cooling_function_data* restrict cooling,
const int convert_part) {
+ /* Initialize the feedback */
+ feedback_init_after_star_formation(sp, e->feedback_props);
+
/* Store the current mass */
const float mass_gas = hydro_get_mass(p);
if (!convert_part) {