diff --git a/.gitignore b/.gitignore
index 5a1118350fb7bf07d2d001b9953bcb963d5d68bd..3c0e7bb7c5582a2809d90534233e752b7eed2681 100644
--- a/.gitignore
+++ b/.gitignore
@@ -62,6 +62,7 @@ examples/Cooling/CoolingRates/cooling_rates
examples/Cooling/CoolingRates/cooling_element_*.dat
examples/Cooling/CoolingRates/cooling_output.dat
examples/SubgridTests/StellarEvolution/StellarEvolutionSolution*
+examples/SubgridTests/CosmologicalStellarEvolution/StellarEvolutionSolution*
tests/testActivePair
tests/testActivePair.sh
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/README b/examples/SubgridTests/CosmologicalStellarEvolution/README
new file mode 100644
index 0000000000000000000000000000000000000000..44a4836e523d4037e5f942e01e444599fa009c1a
--- /dev/null
+++ b/examples/SubgridTests/CosmologicalStellarEvolution/README
@@ -0,0 +1,7 @@
+Example for testing EAGLE stellar feedback. This consists of a uniform box of gas with a star in the center. The amount of feedback can then be checked by summing over the gas particles in the whole box and comparing to the expected amount of feedback from a single star over a given time period.
+
+If only mass enrichment is of interest, the box can be run with ICs generated from a smaller glass (eg glassCube_32.hdf5). In this case, however it is necessary to turn off energy feedback (eg. by setting the return value of compute_SNe in src/stars/EAGLE/stars.h to zero) or using a larger glass (glassCube_64.hdf5).
+
+Use the python script, plot_box_evolution.py to compare total mass evolution of gas particles in the whole box with what is expected based on EAGLE standalone feedback test.
+
+Use plot_paricle_evolution.py to plot the evolution of particles starting in the viscinity of the star at the beginning of the simulation.
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/check_continuous_heating.py b/examples/SubgridTests/CosmologicalStellarEvolution/check_continuous_heating.py
new file mode 100644
index 0000000000000000000000000000000000000000..b5940eba43c89b8af4a883cc8f7022e33293b869
--- /dev/null
+++ b/examples/SubgridTests/CosmologicalStellarEvolution/check_continuous_heating.py
@@ -0,0 +1,137 @@
+# Script for plotting energy evolution of uniform box of gas with single star in the
+# centre when running with stochastic energy injection. It also checks that the change
+# in total energy of the gas particles is within a specified tolerance from what is
+# expected based on the mass of the star particle (Note that this tolerance could be
+# somewhat high because of Poisson noise and the relatively small number of injection
+# events)
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+import os.path
+import numpy as np
+import glob
+
+# 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
+n_elements = 9
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 9,
+'legend.fontsize': 9,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (3.15,3.15),
+'figure.subplot.left' : 0.3,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.18,
+'figure.subplot.top' : 0.92,
+'figure.subplot.wspace' : 0.21,
+'figure.subplot.hspace' : 0.19,
+'lines.markersize' : 6,
+'lines.linewidth' : 2.,
+'text.latex.unicode': True
+}
+
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+# Read the simulation data
+sim = h5py.File("stellar_evolution_0000.hdf5", "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"][0]
+kernel = sim["/HydroScheme"].attrs["Kernel function"][0]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"][0]
+eta = sim["/HydroScheme"].attrs["Kernel eta"][0]
+alpha = sim["/HydroScheme"].attrs["Alpha viscosity"][0]
+H_mass_fraction = sim["/HydroScheme"].attrs["Hydrogen mass fraction"][0]
+H_transition_temp = sim["/HydroScheme"].attrs["Hydrogen ionization transition temperature"][0]
+T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
+T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
+git = sim["Code"].attrs["Git Revision"]
+star_initial_mass = sim["/PartType4/Masses"][0]
+
+# Cosmological parameters
+H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
+gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
+
+# Units
+unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
+unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
+unit_time_in_cgs = sim["/Units"].attrs["Unit time 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_length_in_si = 0.01 * unit_length_in_cgs
+unit_mass_in_si = 0.001 * unit_mass_in_cgs
+unit_time_in_si = unit_time_in_cgs
+
+# Calculate solar mass in internal units
+const_solar_mass = 1.98848e33 / unit_mass_in_cgs
+
+# Define Gyr
+Gyr_in_cgs = 1e9 * 365 * 24 * 3600.
+
+# Find out how many particles (gas and star) we have
+n_parts = sim["/Header"].attrs["NumPart_Total"][0]
+n_sparts = sim["/Header"].attrs["NumPart_Total"][4]
+
+# Declare arrays for data
+masses = zeros((n_parts,n_snapshots))
+star_masses = zeros((n_sparts,n_snapshots))
+internal_energy = zeros((n_parts,n_snapshots))
+velocity_parts = zeros((n_parts,3,n_snapshots))
+time = zeros(n_snapshots)
+
+# Read fields we are checking from snapshots
+#for i in [0,n_snapshots-1]:
+for i in range(n_snapshots):
+ sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
+ print('reading snapshot '+str(i))
+ masses[:,i] = sim["/PartType0/Masses"]
+ internal_energy[:,i] = sim["/PartType0/InternalEnergies"]
+ velocity_parts[:,:,i] = sim["/PartType0/Velocities"]
+ time[i] = sim["/Header"].attrs["Time"][0]
+
+# Check that the total amount of enrichment is as expected.
+# Define tolerance. Note, relatively high value used due to
+# Poisson noise associated with stochastic energy injection.
+eps = 0.15
+
+# Stochastic heating
+vel2 = zeros((n_parts,n_snapshots))
+vel2[:,:] = velocity_parts[:,0,:]*velocity_parts[:,0,:] + velocity_parts[:,1,:]*velocity_parts[:,1,:] + velocity_parts[:,2,:]*velocity_parts[:,2,:]
+total_kinetic_energy_cgs = np.sum(np.multiply(vel2,masses)*0.5,axis = 0) * unit_energy_in_cgs
+total_energy_cgs = np.sum(np.multiply(internal_energy,masses),axis = 0) * unit_energy_in_cgs
+total_energy_released_cgs = total_energy_cgs[n_snapshots-1] - total_energy_cgs[0] + total_kinetic_energy_cgs[n_snapshots-1] - total_kinetic_energy_cgs[0]
+
+# Calculate energy released
+energy_per_sn = 1.0e51 / unit_energy_in_cgs
+SNIa_efficiency = 2.e-3
+SNIa_timescale_Gyr = 2.0
+expected_energy_released_cgs = np.zeros(n_snapshots)
+for i in range(n_snapshots):
+ age_Gyr = time[i] * unit_time_in_cgs / Gyr_in_cgs
+ total_sn = SNIa_efficiency * (1.0 - np.exp(-age_Gyr/SNIa_timescale_Gyr)) / const_solar_mass
+ expected_energy_released_cgs[i] = total_sn * energy_per_sn * unit_energy_in_cgs
+
+# Did we get it right?
+if abs(total_energy_released_cgs - expected_energy_released_cgs[n_snapshots-1])/expected_energy_released_cgs[n_snapshots-1] < eps:
+ print("total stochastic energy release consistent with expectation. total stochastic energy release "+str(total_energy_released_cgs)+" expected "+ str(expected_energy_released_cgs[n_snapshots-1]) + " initial total internal energy "+ str(total_energy_cgs[0] + total_kinetic_energy_cgs[0]))
+else:
+ print("total stochastic energy release "+str(total_energy_released_cgs)+" expected "+ str(expected_energy_released_cgs[n_snapshots-1]) + " initial total internal energy "+ str(total_energy_cgs[0] + total_kinetic_energy_cgs[0]) + " energy change fraction of total " + str(total_energy_released_cgs/(total_energy_cgs[0]+total_kinetic_energy_cgs[0])))
+
+# Plot the energy evolution
+figure()
+subplot(111)
+plot(time*unit_time_in_cgs/Gyr_in_cgs, total_energy_cgs + total_kinetic_energy_cgs - total_energy_cgs[0] - total_kinetic_energy_cgs[0],color='k', linewidth=0.5, label="SWIFT")
+plot(time*unit_time_in_cgs/Gyr_in_cgs, expected_energy_released_cgs,color = 'r', linewidth=0.5, label="expected")
+xlabel("Time (Gyr)")
+ylabel("Total energy (erg)")
+legend()
+savefig("continuous_energy_evolution.png", dpi=200)
+
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/check_stellar_evolution.py b/examples/SubgridTests/CosmologicalStellarEvolution/check_stellar_evolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..5680eb4d64f29ab32831d995e31ae9c27de82a71
--- /dev/null
+++ b/examples/SubgridTests/CosmologicalStellarEvolution/check_stellar_evolution.py
@@ -0,0 +1,318 @@
+import matplotlib
+
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+import os.path
+import numpy as np
+import glob
+
+# 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
+)
+n_elements = 9
+
+# Plot parameters
+params = {
+ "axes.labelsize": 10,
+ "axes.titlesize": 10,
+ "font.size": 9,
+ "legend.fontsize": 9,
+ "xtick.labelsize": 10,
+ "ytick.labelsize": 10,
+ "text.usetex": True,
+ "figure.figsize": (3.15, 3.15),
+ "figure.subplot.left": 0.3,
+ "figure.subplot.right": 0.99,
+ "figure.subplot.bottom": 0.18,
+ "figure.subplot.top": 0.92,
+ "figure.subplot.wspace": 0.21,
+ "figure.subplot.hspace": 0.19,
+ "lines.markersize": 6,
+ "lines.linewidth": 2.0,
+ "text.latex.unicode": True,
+}
+
+rcParams.update(params)
+rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
+
+# Read the simulation data
+sim = h5py.File("stellar_evolution_0000.hdf5", "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"][0]
+kernel = sim["/HydroScheme"].attrs["Kernel function"][0]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"][0]
+eta = sim["/HydroScheme"].attrs["Kernel eta"][0]
+alpha = sim["/HydroScheme"].attrs["Alpha viscosity"][0]
+H_mass_fraction = sim["/HydroScheme"].attrs["Hydrogen mass fraction"][0]
+H_transition_temp = sim["/HydroScheme"].attrs[
+ "Hydrogen ionization transition temperature"
+][0]
+T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
+T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
+git = sim["Code"].attrs["Git Revision"]
+star_initial_mass = sim["/PartType4/Masses"][0]
+
+# Cosmological parameters
+H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
+gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
+
+# Units
+unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
+unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
+unit_time_in_cgs = sim["/Units"].attrs["Unit time 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_length_in_si = 0.01 * unit_length_in_cgs
+unit_mass_in_si = 0.001 * unit_mass_in_cgs
+unit_time_in_si = unit_time_in_cgs
+
+# Find out how many particles (gas and star) we have
+n_parts = sim["/Header"].attrs["NumPart_Total"][0]
+n_sparts = sim["/Header"].attrs["NumPart_Total"][4]
+
+# Declare arrays for data
+masses = zeros((n_parts, n_snapshots))
+star_masses = zeros((n_sparts, n_snapshots))
+mass_from_AGB = zeros((n_parts, n_snapshots))
+metal_mass_frac_from_AGB = zeros((n_parts, n_snapshots))
+mass_from_SNII = zeros((n_parts, n_snapshots))
+metal_mass_frac_from_SNII = zeros((n_parts, n_snapshots))
+mass_from_SNIa = zeros((n_parts, n_snapshots))
+metal_mass_frac_from_SNIa = zeros((n_parts, n_snapshots))
+iron_mass_frac_from_SNIa = zeros((n_parts, n_snapshots))
+metallicity = zeros((n_parts, n_snapshots))
+abundances = zeros((n_parts, n_elements, n_snapshots))
+internal_energy = zeros((n_parts, n_snapshots))
+coord_parts = zeros((n_parts, 3))
+velocity_parts = zeros((n_parts, 3, n_snapshots))
+coord_sparts = zeros(3)
+time = zeros(n_snapshots)
+
+# Read fields we are checking from snapshots
+for i in range(n_snapshots):
+ sim = h5py.File("stellar_evolution_%04d.hdf5" % i, "r")
+ print("reading snapshot " + str(i))
+ abundances[:, :, i] = sim["/PartType0/ElementMassFractions"]
+ metallicity[:, i] = sim["/PartType0/Metallicities"]
+ masses[:, i] = sim["/PartType0/Masses"]
+ star_masses[:, i] = sim["/PartType4/Masses"]
+ mass_from_AGB[:, i] = sim["/PartType0/TotalMassFromAGB"]
+ metal_mass_frac_from_AGB[:, i] = sim["/PartType0/MetalMassFracFromAGB"]
+ mass_from_SNII[:, i] = sim["/PartType0/TotalMassFromSNII"]
+ metal_mass_frac_from_SNII[:, i] = sim["/PartType0/MetalMassFracFromSNII"]
+ mass_from_SNIa[:, i] = sim["/PartType0/TotalMassFromSNIa"]
+ metal_mass_frac_from_SNIa[:, i] = sim["/PartType0/MetalMassFracFromSNIa"]
+ iron_mass_frac_from_SNIa[:, i] = sim["/PartType0/IronMassFracFromSNIa"]
+ internal_energy[:, i] = sim["/PartType0/InternalEnergies"]
+ velocity_parts[:, :, i] = sim["/PartType0/Velocities"]
+ time[i] = sim["/Header"].attrs["Time"][0]
+
+# Define ejecta factor
+ejecta_factor = 1.0e-2
+ejecta_factor_metallicity = 1.0 - 2.0 / n_elements
+ejecta_factor_abundances = 1.0 / n_elements
+ejected_mass = star_initial_mass
+energy_per_SNe = 1.0e51 / unit_energy_in_cgs
+
+# Check that the total amount of enrichment is as expected.
+# Define tolerance
+eps = 0.01
+
+# Total mass
+total_part_mass = np.sum(masses, axis=0)
+if (
+ abs(
+ (total_part_mass[n_snapshots - 1] - total_part_mass[0]) / total_part_mass[0]
+ - ejected_mass / total_part_mass[0]
+ )
+ * total_part_mass[0]
+ / ejected_mass
+ < eps
+):
+ print("total mass released consistent with expectation")
+else:
+ print(
+ "mass increase "
+ + str(total_part_mass[n_snapshots - 1] / total_part_mass[0])
+ + " expected "
+ + str(1.0 + ejected_mass / total_part_mass[0])
+ )
+
+# Check that mass is conserved (i.e. total star mass decreases by same amount as total gas mass increases)
+total_spart_mass = np.sum(star_masses, axis=0)
+if (
+ abs(
+ (total_part_mass[n_snapshots - 1] + total_spart_mass[n_snapshots - 1])
+ / (total_part_mass[0] + total_spart_mass[0])
+ - 1.0
+ )
+ < eps ** 3
+):
+ print("total mass conserved")
+else:
+ print(
+ "initial part, spart mass "
+ + str(total_part_mass[0])
+ + " "
+ + str(total_spart_mass[0])
+ + " final mass "
+ + str(total_part_mass[n_snapshots - 1])
+ + " "
+ + str(total_spart_mass[n_snapshots - 1])
+ )
+
+# Total metal mass from AGB
+total_metal_mass_AGB = np.sum(np.multiply(metal_mass_frac_from_AGB, masses), axis=0)
+expected_metal_mass_AGB = ejecta_factor * ejected_mass
+if (
+ abs(total_metal_mass_AGB[n_snapshots - 1] - expected_metal_mass_AGB)
+ / expected_metal_mass_AGB
+ < eps
+):
+ print("total AGB metal mass released consistent with expectation")
+else:
+ print(
+ "total AGB metal mass "
+ + str(total_metal_mass_AGB[n_snapshots - 1])
+ + " expected "
+ + str(expected_metal_mass_AGB)
+ )
+
+# Total mass from AGB
+total_AGB_mass = np.sum(mass_from_AGB, axis=0)
+expected_AGB_mass = ejecta_factor * ejected_mass
+if abs(total_AGB_mass[n_snapshots - 1] - expected_AGB_mass) / expected_AGB_mass < eps:
+ print("total AGB mass released consistent with expectation")
+else:
+ print(
+ "total AGB mass "
+ + str(total_AGB_mass[n_snapshots - 1])
+ + " expected "
+ + str(expected_AGB_mass)
+ )
+
+# Total metal mass from SNII
+total_metal_mass_SNII = np.sum(np.multiply(metal_mass_frac_from_SNII, masses), axis=0)
+expected_metal_mass_SNII = ejecta_factor * ejected_mass
+if (
+ abs(total_metal_mass_SNII[n_snapshots - 1] - expected_metal_mass_SNII)
+ / expected_metal_mass_SNII
+ < eps
+):
+ print("total SNII metal mass released consistent with expectation")
+else:
+ print(
+ "total SNII metal mass "
+ + str(total_metal_mass_SNII[n_snapshots - 1])
+ + " expected "
+ + str(expected_metal_mass_SNII)
+ )
+
+# Total mass from SNII
+total_SNII_mass = np.sum(mass_from_SNII, axis=0)
+expected_SNII_mass = ejecta_factor * ejected_mass
+if (
+ abs(total_SNII_mass[n_snapshots - 1] - expected_SNII_mass) / expected_SNII_mass
+ < eps
+):
+ print("total SNII mass released consistent with expectation")
+else:
+ print(
+ "total SNII mass "
+ + str(total_SNII_mass[n_snapshots - 1])
+ + " expected "
+ + str(expected_SNII_mass)
+ )
+
+# Total metal mass from SNIa
+total_metal_mass_SNIa = np.sum(np.multiply(metal_mass_frac_from_SNIa, masses), axis=0)
+expected_metal_mass_SNIa = ejecta_factor * ejected_mass
+if (
+ abs(total_metal_mass_SNIa[n_snapshots - 1] - expected_metal_mass_SNIa)
+ / expected_metal_mass_SNIa
+ < eps
+):
+ print("total SNIa metal mass released consistent with expectation")
+else:
+ print(
+ "total SNIa metal mass "
+ + str(total_metal_mass_SNIa[n_snapshots - 1])
+ + " expected "
+ + str(expected_metal_mass_SNIa)
+ )
+
+# Total iron mass from SNIa
+total_iron_mass_SNIa = np.sum(np.multiply(iron_mass_frac_from_SNIa, masses), axis=0)
+expected_iron_mass_SNIa = ejecta_factor * ejected_mass
+if (
+ abs(total_iron_mass_SNIa[n_snapshots - 1] - expected_iron_mass_SNIa)
+ / expected_iron_mass_SNIa
+ < eps
+):
+ print("total SNIa iron mass released consistent with expectation")
+else:
+ print(
+ "total SNIa iron mass "
+ + str(total_iron_mass_SNIa[n_snapshots - 1])
+ + " expected "
+ + str(expected_iron_mass_SNIa)
+ )
+
+# Total mass from SNIa
+total_SNIa_mass = np.sum(mass_from_SNIa, axis=0)
+expected_SNIa_mass = ejecta_factor * ejected_mass
+if (
+ abs(total_SNIa_mass[n_snapshots - 1] - expected_SNIa_mass) / expected_SNIa_mass
+ < eps
+):
+ print("total SNIa mass released consistent with expectation")
+else:
+ print(
+ "total SNIa mass "
+ + str(total_SNIa_mass[n_snapshots - 1])
+ + " expected "
+ + str(expected_SNIa_mass)
+ )
+
+# Total metal mass
+total_metal_mass = np.sum(np.multiply(metallicity, masses), axis=0)
+expected_metal_mass = ejecta_factor_metallicity * ejected_mass
+if (
+ abs(total_metal_mass[n_snapshots - 1] - expected_metal_mass) / expected_metal_mass
+ < eps
+):
+ print("total metal mass released consistent with expectation")
+else:
+ print(
+ "total metal mass "
+ + str(total_metal_mass[n_snapshots - 1])
+ + " expected "
+ + str(expected_metal_mass)
+ )
+
+# Total mass for each element
+expected_element_mass = ejecta_factor_abundances * ejected_mass
+for i in range(n_elements):
+ total_element_mass = np.sum(np.multiply(abundances[:, i, :], masses), axis=0)
+ if (
+ abs(total_element_mass[n_snapshots - 1] - expected_element_mass)
+ / expected_element_mass
+ < eps
+ ):
+ print(
+ "total element mass released consistent with expectation for element "
+ + str(i)
+ )
+ else:
+ print(
+ "total element mass "
+ + str(total_element_mass[n_snapshots - 1])
+ + " expected "
+ + str(expected_element_mass)
+ + " for element "
+ + str(i)
+ )
+
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/check_stochastic_heating.py b/examples/SubgridTests/CosmologicalStellarEvolution/check_stochastic_heating.py
new file mode 100644
index 0000000000000000000000000000000000000000..1cacc13653d821da3abd2a09566be347608c64f7
--- /dev/null
+++ b/examples/SubgridTests/CosmologicalStellarEvolution/check_stochastic_heating.py
@@ -0,0 +1,137 @@
+# Script for plotting energy evolution of uniform box of gas with single star in the
+# centre when running with stochastic energy injection. It also checks that the change
+# in total energy of the gas particles is within a specified tolerance from what is
+# expected based on the mass of the star particle (Note that this tolerance could be
+# somewhat high because of Poisson noise and the relatively small number of injection
+# events)
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+import os.path
+import numpy as np
+import glob
+
+# 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
+n_elements = 9
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 9,
+'legend.fontsize': 9,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (3.15,3.15),
+'figure.subplot.left' : 0.3,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.18,
+'figure.subplot.top' : 0.92,
+'figure.subplot.wspace' : 0.21,
+'figure.subplot.hspace' : 0.19,
+'lines.markersize' : 6,
+'lines.linewidth' : 2.,
+'text.latex.unicode': True
+}
+
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+# Read the simulation data
+sim = h5py.File("stellar_evolution_0000.hdf5", "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"][0]
+kernel = sim["/HydroScheme"].attrs["Kernel function"][0]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"][0]
+eta = sim["/HydroScheme"].attrs["Kernel eta"][0]
+alpha = sim["/HydroScheme"].attrs["Alpha viscosity"][0]
+H_mass_fraction = sim["/HydroScheme"].attrs["Hydrogen mass fraction"][0]
+H_transition_temp = sim["/HydroScheme"].attrs["Hydrogen ionization transition temperature"][0]
+T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
+T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
+git = sim["Code"].attrs["Git Revision"]
+star_initial_mass = sim["/PartType4/Masses"][0]
+
+# Cosmological parameters
+H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
+gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
+
+# Units
+unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
+unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
+unit_time_in_cgs = sim["/Units"].attrs["Unit time 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_length_in_si = 0.01 * unit_length_in_cgs
+unit_mass_in_si = 0.001 * unit_mass_in_cgs
+unit_time_in_si = unit_time_in_cgs
+
+# Calculate solar mass in internal units
+const_solar_mass = 1.98848e33 / unit_mass_in_cgs
+
+# Define Gyr
+Gyr_in_cgs = 1e9 * 365 * 24 * 3600.
+
+# Find out how many particles (gas and star) we have
+n_parts = sim["/Header"].attrs["NumPart_Total"][0]
+n_sparts = sim["/Header"].attrs["NumPart_Total"][4]
+
+# Declare arrays for data
+masses = zeros((n_parts,n_snapshots))
+star_masses = zeros((n_sparts,n_snapshots))
+internal_energy = zeros((n_parts,n_snapshots))
+velocity_parts = zeros((n_parts,3,n_snapshots))
+time = zeros(n_snapshots)
+
+# Read fields we are checking from snapshots
+#for i in [0,n_snapshots-1]:
+for i in range(n_snapshots):
+ sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
+ print('reading snapshot '+str(i))
+ masses[:,i] = sim["/PartType0/Masses"]
+ internal_energy[:,i] = sim["/PartType0/InternalEnergies"]
+ velocity_parts[:,:,i] = sim["/PartType0/Velocities"]
+ time[i] = sim["/Header"].attrs["Time"][0]
+
+# Check that the total amount of enrichment is as expected.
+# Define tolerance. Note, relatively high value used due to
+# Poisson noise associated with stochastic energy injection.
+eps = 0.15
+
+# Stochastic heating
+vel2 = zeros((n_parts,n_snapshots))
+vel2[:,:] = velocity_parts[:,0,:]*velocity_parts[:,0,:] + velocity_parts[:,1,:]*velocity_parts[:,1,:] + velocity_parts[:,2,:]*velocity_parts[:,2,:]
+total_kinetic_energy_cgs = np.sum(np.multiply(vel2,masses)*0.5,axis = 0) * unit_energy_in_cgs
+total_energy_cgs = np.sum(np.multiply(internal_energy,masses),axis = 0) * unit_energy_in_cgs
+total_energy_released_cgs = total_energy_cgs[n_snapshots-1] - total_energy_cgs[0] + total_kinetic_energy_cgs[n_snapshots-1] - total_kinetic_energy_cgs[0]
+
+# Calculate energy released
+num_SNII_per_msun = 1.73621e-02
+energy_per_sn = 1.0e51 / unit_energy_in_cgs
+expected_energy_released_cgs = np.zeros(n_snapshots)
+for i in range(n_snapshots):
+ if time[i]*unit_time_in_cgs/Gyr_in_cgs < 0.03:
+ expected_energy_released_cgs[i] = 0
+ else:
+ expected_energy_released_cgs[i] = num_SNII_per_msun * star_initial_mass / const_solar_mass * energy_per_sn * unit_energy_in_cgs
+
+# Did we get it right?
+if abs(total_energy_released_cgs - expected_energy_released_cgs[n_snapshots-1])/expected_energy_released_cgs[n_snapshots-1] < eps:
+ print("total stochastic energy release consistent with expectation. total stochastic energy release "+str(total_energy_released_cgs)+" expected "+ str(expected_energy_released_cgs[n_snapshots-1]) + " initial total internal energy "+ str(total_energy_cgs[0] + total_kinetic_energy_cgs[0]))
+else:
+ print("total stochastic energy release "+str(total_energy_released_cgs)+" expected "+ str(expected_energy_released_cgs[n_snapshots-1]) + " initial total internal energy "+ str(total_energy_cgs[0] + total_kinetic_energy_cgs[0]) + " energy change fraction of total " + str(total_energy_released_cgs/(total_energy_cgs[0]+total_kinetic_energy_cgs[0])))
+
+# Plot the energy evolution
+figure()
+subplot(111)
+plot(time*unit_time_in_cgs/Gyr_in_cgs, total_energy_cgs + total_kinetic_energy_cgs - total_energy_cgs[0] - total_kinetic_energy_cgs[0],color='k', linewidth=0.5, label="SWIFT")
+plot(time*unit_time_in_cgs/Gyr_in_cgs, expected_energy_released_cgs,color = 'r', linewidth=0.5, label="expected")
+xlabel("Time (Gyr)")
+ylabel("Total energy (erg)")
+legend()
+savefig("stochastic_energy_evolution.png", dpi=200)
+
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/getEagleYieldTable.sh b/examples/SubgridTests/CosmologicalStellarEvolution/getEagleYieldTable.sh
new file mode 100755
index 0000000000000000000000000000000000000000..26eef020cab82acee2c80e88089df1790b281eab
--- /dev/null
+++ b/examples/SubgridTests/CosmologicalStellarEvolution/getEagleYieldTable.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/YieldTables/EAGLE/yieldtables.tar.gz
+tar -xf yieldtables.tar.gz
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/getGlass.sh b/examples/SubgridTests/CosmologicalStellarEvolution/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..ffd92e88deae6e91237059adac2a6c2067caee46
--- /dev/null
+++ b/examples/SubgridTests/CosmologicalStellarEvolution/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/getSolutions.sh b/examples/SubgridTests/CosmologicalStellarEvolution/getSolutions.sh
new file mode 100755
index 0000000000000000000000000000000000000000..1fe0f1507a7efa1da843970ddcede681e846e4ec
--- /dev/null
+++ b/examples/SubgridTests/CosmologicalStellarEvolution/getSolutions.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/StellarEvolutionSolution.tar.gz
+tar -xvzf StellarEvolutionSolution.tar.gz
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/makeIC.py b/examples/SubgridTests/CosmologicalStellarEvolution/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..abc2d8126ae8ac1173f8918ffc818f9cb6bcb4fe
--- /dev/null
+++ b/examples/SubgridTests/CosmologicalStellarEvolution/makeIC.py
@@ -0,0 +1,144 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+import h5py
+from numpy import *
+
+# Some constants
+solar_mass_cgs = 1.988480e33
+kpc_in_cm = 3.085678e21
+mp_cgs = 1.67e-24
+boltzmann_k_cgs = 1.38e-16
+
+# Parameters
+gamma = 5./3. # Gas adiabatic index
+rho_cgs = mp_cgs # Background density
+u0_cgs = 1.2e12 # Desired initial internal energy (1.2e12 ~ 10^4K)
+P_cgs = rho_cgs*u0_cgs*(gamma - 1.) # Background pressure
+fileName = "stellar_evolution.hdf5"
+
+# Units
+unit_l_cgs = 3.085678e24 # kpc
+unit_m_cgs = 1.988480e43 # 10^10 Msun
+unit_v_cgs = 1e5 # km / s
+unit_A_cgs = 1.
+unit_T_cgs = 1.
+unit_t_cgs = unit_l_cgs / unit_v_cgs
+
+boxsize_cgs = 10. * kpc_in_cm
+vol_cgs = boxsize_cgs**3
+
+#---------------------------------------------------
+glass = h5py.File("glassCube_32.hdf5", "r")
+
+# Read particle positions and h from the glass
+pos = glass["/PartType0/Coordinates"][:,:]
+eps = 1e-6
+pos = (pos - pos.min()) / (pos.max() - pos.min() + eps) * boxsize_cgs / unit_l_cgs
+h = glass["/PartType0/SmoothingLength"][:] * boxsize_cgs / unit_l_cgs
+
+numPart = size(h)
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m_cgs = zeros(numPart)
+u_cgs = zeros(numPart)
+m = zeros(numPart)
+u = zeros(numPart)
+
+m_cgs[:] = rho_cgs * vol_cgs / numPart
+u_cgs[:] = P_cgs / (rho_cgs * (gamma - 1))
+
+# Stars
+star_pos = zeros((1, 3))
+star_pos[:,:] = 0.5 * boxsize_cgs / unit_l_cgs
+
+star_v = zeros((1, 3))
+star_v[:,:] = 0.
+
+# increase mass to keep it at center
+star_m_cgs = m_cgs[0]
+star_ids = array([numPart + 1])
+star_h = array([h.max()])
+
+#--------------------------------------------------
+
+# Check quantities are correct for debugging
+print("part mass/msun " + str(m_cgs[0]/solar_mass_cgs) + " stellar mass/msun " + str(star_m_cgs/solar_mass_cgs))
+print("boxsize kpc " + str(boxsize_cgs/kpc_in_cm))
+print("density cm^-3 " + str(rho_cgs/mp_cgs))
+print("initial temperature K " + str(u_cgs[0] / boltzmann_k_cgs*((gamma - 1)*rho_cgs)))
+
+# Convert to internal units
+star_m = star_m_cgs/unit_m_cgs
+m[:] = m_cgs/unit_m_cgs
+u[:] = u_cgs*unit_v_cgs**-2
+boxsize = boxsize_cgs/unit_l_cgs
+
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [boxsize]*3
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 1, 0]
+#grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 1, 0]
+#grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+grp.attrs["Dimension"] = 3
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 0
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = unit_l_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = unit_m_cgs
+grp.attrs["Unit time in cgs (U_t)"] = unit_t_cgs
+grp.attrs["Unit current in cgs (U_I)"] = unit_A_cgs
+grp.attrs["Unit temperature in cgs (U_T)"] = unit_T_cgs
+
+#Particle group
+grp = file.create_group("/PartType0")
+grp.create_dataset('Coordinates', data=pos, dtype='d')
+grp.create_dataset('Velocities', data=v, dtype='f')
+grp.create_dataset('Masses', data=m, dtype='f')
+grp.create_dataset('SmoothingLength', data=h, dtype='f')
+grp.create_dataset('InternalEnergy', data=u, dtype='f')
+grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+
+# stellar group
+grp = file.create_group("/PartType4")
+grp.create_dataset("Coordinates", data=star_pos, dtype="d")
+grp.create_dataset('Velocities', data=star_v, dtype='f')
+grp.create_dataset('Masses', data=star_m, dtype='f')
+grp.create_dataset('SmoothingLength', data=star_h, dtype='f')
+grp.create_dataset('ParticleIDs', data=star_ids, dtype='L')
+
+file.close()
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/plot_box_evolution.py b/examples/SubgridTests/CosmologicalStellarEvolution/plot_box_evolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..96f4acf4352660d9ff05ffa62b83b156e654da9f
--- /dev/null
+++ b/examples/SubgridTests/CosmologicalStellarEvolution/plot_box_evolution.py
@@ -0,0 +1,253 @@
+###############################################################################
+ # 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
+import h5py
+import numpy as np
+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
+}
+rcParams.update(params)
+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
+n_elements = 9
+
+# Read the simulation data
+sim = h5py.File("stellar_evolution_0000.hdf5", "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+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
+
+# Units
+unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
+unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
+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.
+Msun_in_cgs = 1.98848e33
+
+# Declare arrays to store SWIFT data
+swift_box_gas_mass = zeros(n_snapshots)
+swift_box_star_mass = zeros(n_snapshots)
+swift_box_gas_metal_mass = zeros(n_snapshots)
+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.
+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]
+
+ masses = sim["/PartType0/Masses"][:]
+ swift_box_gas_mass[i] = np.sum(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)
+
+ 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]
+
+ if i == 0:
+ swift_mean_u_start = np.mean(u)
+
+ 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
+
+# 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
+
+# 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
+
+
+
+# Plot the interesting quantities
+figure()
+
+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)
+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)
+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']
+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)
+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)
+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")
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/plot_particle_evolution.py b/examples/SubgridTests/CosmologicalStellarEvolution/plot_particle_evolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..be1588f9b707d448b2905611defd9e760c5f91de
--- /dev/null
+++ b/examples/SubgridTests/CosmologicalStellarEvolution/plot_particle_evolution.py
@@ -0,0 +1,223 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+# 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/>.
+#
+##############################################################################
+
+# Assuming output snapshots contain evolution of box of gas with star at its
+# centre, this script will plot the evolution of the radial velocities, internal
+# energies, mass and metallicities of the nearest n particles to the star over
+# the duration of the simulation.
+
+import matplotlib
+
+matplotlib.use("Agg")
+from pylab import *
+from scipy import stats
+import h5py
+import numpy as np
+import glob
+import os.path
+
+# Function to find index in array a for each element in array b
+def find_indices(a, b):
+ result = np.zeros(len(b))
+ for i in range(len(b)):
+ result[i] = ((np.where(a == b[i]))[0])[0]
+ return result
+
+
+# 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.1,
+ "figure.subplot.right": 0.99,
+ "figure.subplot.bottom": 0.1,
+ "figure.subplot.top": 0.95,
+ "figure.subplot.wspace": 0.2,
+ "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"]})
+
+
+# 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
+)
+n_particles_to_plot = 500
+
+# Read the simulation data
+sim = h5py.File("stellar_evolution_0000.hdf5", "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = sim["Code"].attrs["Git Revision"]
+
+# Units
+unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
+unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
+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_length_in_si = 0.01 * unit_length_in_cgs
+unit_mass_in_si = 0.001 * unit_mass_in_cgs
+unit_time_in_si = unit_time_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
+Myr_in_cgs = 3.154e13
+Msun_in_cgs = 1.989e33
+
+# Read data of zeroth snapshot
+pos = sim["/PartType0/Coordinates"][:, :]
+x = pos[:, 0] - boxSize / 2
+y = pos[:, 1] - boxSize / 2
+z = pos[:, 2] - boxSize / 2
+vel = sim["/PartType0/Velocities"][:, :]
+r = sqrt(x ** 2 + y ** 2 + z ** 2)
+v_r = (x * vel[:, 0] + y * vel[:, 1] + z * vel[:, 2]) / r
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
+mass = sim["/PartType0/Masses"][:]
+IDs = sim["/PartType0/ParticleIDs"][:]
+
+# Find which particles are closest to centre of box
+index = argsort(r)
+part_IDs_to_plot = zeros(n_particles_to_plot)
+part_IDs_to_plot = np.sort(IDs[index[0:n_particles_to_plot]])
+
+# Declare arrrays to plot
+masses_to_plot = zeros((n_particles_to_plot, n_snapshots))
+v_r_to_plot = zeros((n_particles_to_plot, n_snapshots))
+metallicities_to_plot = zeros((n_particles_to_plot, n_snapshots))
+internal_energies_to_plot = zeros((n_particles_to_plot, n_snapshots))
+t = zeros(n_snapshots)
+
+# Read data from rest of snapshots
+for i in range(n_snapshots):
+ print("reading snapshot " + str(i))
+ # Read the simulation data
+ sim = h5py.File("stellar_evolution_%04d.hdf5" % i, "r")
+ t[i] = sim["/Header"].attrs["Time"][0]
+
+ pos = sim["/PartType0/Coordinates"][:, :]
+ x = pos[:, 0] - boxSize / 2
+ y = pos[:, 1] - boxSize / 2
+ z = pos[:, 2] - boxSize / 2
+ vel = sim["/PartType0/Velocities"][:, :]
+ r = sqrt(x ** 2 + y ** 2 + z ** 2)
+ v_r = (x * vel[:, 0] + y * vel[:, 1] + z * vel[:, 2]) / r
+ u = sim["/PartType0/InternalEnergies"][:]
+ S = sim["/PartType0/Entropies"][:]
+ P = sim["/PartType0/Pressures"][:]
+ rho = sim["/PartType0/Densities"][:]
+ mass = sim["/PartType0/Masses"][:]
+ metallicity = sim["/PartType0/Metallicities"][:]
+ internal_energy = sim["/PartType0/InternalEnergies"][:]
+ IDs = sim["/PartType0/ParticleIDs"][:]
+
+ # Find which particles we want to plot and store their data
+ indices = (find_indices(IDs, part_IDs_to_plot)).astype(int)
+ masses_to_plot[:, i] = mass[indices[:]]
+ v_r_to_plot[:, i] = v_r[indices[:]]
+ metallicities_to_plot[:, i] = metallicity[indices[:]]
+ internal_energies_to_plot[:, i] = internal_energy[indices[:]]
+
+
+# Plot the interesting quantities
+figure()
+
+# Radial velocity --------------------------------
+subplot(221)
+for j in range(n_particles_to_plot):
+ plot(
+ t * unit_time_in_cgs / Myr_in_cgs,
+ v_r_to_plot[j, :] * unit_vel_in_cgs,
+ linewidth=0.5,
+ color="k",
+ ms=0.5,
+ alpha=0.1,
+ )
+xlabel("Time (Myr)", labelpad=0)
+ylabel("Radial velocity $(\\rm{cm} \cdot \\rm{s}^{-1})$", labelpad=0)
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
+
+# Internal energy --------------------------------
+subplot(222)
+for j in range(n_particles_to_plot):
+ plot(
+ t * unit_time_in_cgs / Myr_in_cgs,
+ internal_energies_to_plot[j, :] * unit_energy_in_cgs / unit_mass_in_cgs,
+ linewidth=0.5,
+ color="k",
+ ms=0.5,
+ alpha=0.1,
+ )
+xlabel("Time (Myr)", labelpad=0)
+ylabel("Internal energy $(\\rm{erg} \cdot \\rm{g}^{-1})$", labelpad=2)
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
+
+# Masses --------------------------------
+subplot(223)
+for j in range(n_particles_to_plot):
+ plot(
+ t * unit_time_in_cgs / Myr_in_cgs,
+ masses_to_plot[j, :] * unit_mass_in_cgs / Msun_in_cgs,
+ linewidth=0.5,
+ color="k",
+ ms=0.5,
+ alpha=0.1,
+ )
+xlabel("Time (Myr)", labelpad=0)
+ylabel("Mass (Msun)", labelpad=2)
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
+
+# Metallicities --------------------------------
+subplot(224)
+for j in range(n_particles_to_plot):
+ plot(
+ t * unit_time_in_cgs / Myr_in_cgs,
+ metallicities_to_plot[j, :],
+ linewidth=0.5,
+ color="k",
+ ms=0.5,
+ alpha=0.1,
+ )
+xlabel("Time (Myr)", labelpad=0)
+ylabel("Metallicity", labelpad=2)
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
+
+savefig("particle_evolution.png", dpi=200)
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/run.sh b/examples/SubgridTests/CosmologicalStellarEvolution/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..637a100567f691d883b95f78c4513e0640c1932a
--- /dev/null
+++ b/examples/SubgridTests/CosmologicalStellarEvolution/run.sh
@@ -0,0 +1,47 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e glassCube_32.hdf5 ]
+then
+ echo "Fetching initial glass file for the Supernovae feedback example..."
+ ./getGlass.sh
+fi
+if [ ! -e stellar_evolution.hdf5 ]
+then
+ echo "Generating initial conditions for the 3D stellar evolution example..."
+ python makeIC.py
+fi
+
+# Get the Yield tables
+if [ ! -e yieldtables ]
+then
+ echo "Fetching Yield tables..."
+ ./getEagleYieldTable.sh
+fi
+
+# Get the solutions
+if [ ! -e StellarEvolutionSolution ]
+then
+ echo "Fetching solutions ..."
+ ./getSolutions.sh
+fi
+
+../../swift --feedback --stars --hydro --cosmology --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 --cosmology --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 --cosmology --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 --cosmology --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 --cosmology --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
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/stellar_evolution.yml b/examples/SubgridTests/CosmologicalStellarEvolution/stellar_evolution.yml
new file mode 100644
index 0000000000000000000000000000000000000000..9b8c3e34dad20eba560a7316f16364b76a088c05
--- /dev/null
+++ b/examples/SubgridTests/CosmologicalStellarEvolution/stellar_evolution.yml
@@ -0,0 +1,120 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
+ UnitLength_in_cgs: 3.08567758e24 # Mpc in centimeters
+ UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
+ UnitCurrent_in_cgs: 1. # Amperes
+ UnitTemp_in_cgs: 1. # Kelvin
+
+# Cosmological parameters
+Cosmology:
+ h: 0.6777 # Reduced Hubble constant
+ a_begin: 0.0099 # Initial scale-factor of the simulation (z = 100.0)
+ a_end: 1.0 # Final scale factor of the simulation
+ Omega_m: 0.307 # Matter density parameter
+ Omega_lambda: 0.693 # Dark-energy density parameter
+ Omega_b: 0.0482519 # Baryon density parameter
+
+# Parameters governing the time integration
+TimeIntegration:
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 5e-3 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: stellar_evolution # Common part of the name of output files
+ delta_time: 1.03
+ scale_factor_first: 0.02
+ compression: 4
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ scale_factor_first: 0.00991
+ delta_time: 1.1
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100. # Kelvin
+
+# Properties of the stars
+Stars:
+ birth_time: 0.00991 # Give the star in the ICs a decent birth time
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./stellar_evolution.hdf5 # The file to read
+ periodic: 1
+
+Scheduler:
+ max_top_level_cells: 8
+
+# Parameters for the EAGLE "equation of state"
+EAGLEEntropyFloor:
+ Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
+ Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
+ Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
+ Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
+ Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
+ Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
+
+# Metallicites read in for the gas and star
+EAGLEChemistry:
+ init_abundance_metal: 0.01
+ init_abundance_Hydrogen: 0.752
+ init_abundance_Helium: 0.248
+ init_abundance_Carbon: 0.0
+ init_abundance_Nitrogen: 0.0
+ init_abundance_Oxygen: 0.0
+ init_abundance_Neon: 0.0
+ init_abundance_Magnesium: 0.0
+ init_abundance_Silicon: 0.0
+ init_abundance_Iron: 0.0
+
+# Standard EAGLE cooling options
+EAGLECooling:
+ dir_name: ./coolingtables/ # Location of the Wiersma+08 cooling tables
+ H_reion_z: 11.5 # Redshift of Hydrogen re-ionization
+ H_reion_eV_p_H: 2.0
+ He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
+ He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
+ He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
+
+# Properties of the EAGLE feedback and enrichment model.
+EAGLEFeedback:
+ use_SNII_feedback: 0 # Global switch for SNII thermal (stochastic) 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.
+ filename: ./yieldtables/ # Path to the directory containing the EAGLE yield tables.
+ IMF_min_mass_Msun: 0.1 # Minimal stellar mass considered for the Chabrier IMF in solar masses.
+ 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 feedback (not SNII enrichment!) in solar masses.
+ SNII_max_mass_Msun: 100.0 # Maximal mass considered for SNII feedback (not SNII enrichment!) in solar masses.
+ SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event.
+ SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
+ SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
+ SNII_energy_fraction_min: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_max: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
+ 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.
+ SNIa_max_mass_Msun: 8.0 # Maximal mass considered for SNIa feedback and enrichment in solar masses.
+ SNIa_timescale_Gyr: 2.0 # Time-scale of the exponential decay of the SNIa rates in Gyr.
+ SNIa_efficiency_p_Msun: 0.002 # Normalisation of the SNIa rates in inverse solar masses.
+ SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
+ AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
+ SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
+ SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
+ SNII_yield_factor_Nitrogen: 1.0 # (Optional) Correction factor to apply to the Nitrogen yield from the SNII channel.
+ SNII_yield_factor_Oxygen: 1.0 # (Optional) Correction factor to apply to the Oxygen yield from the SNII channel.
+ SNII_yield_factor_Neon: 1.0 # (Optional) Correction factor to apply to the Neon yield from the SNII channel.
+ SNII_yield_factor_Magnesium: 2.0 # (Optional) Correction factor to apply to the Magnesium yield from the SNII channel.
+ SNII_yield_factor_Silicon: 1.0 # (Optional) Correction factor to apply to the Silicon yield from the SNII channel.
+ SNII_yield_factor_Iron: 0.5 # (Optional) Correction factor to apply to the Iron yield from the SNII channel.
diff --git a/examples/SubgridTests/StellarEvolution/run.sh b/examples/SubgridTests/StellarEvolution/run.sh
index 7cf174cbfca12dad84eca72ab0329eb065f1c046..ec427c4ead207749041002d46670bfc91a4e4d78 100755
--- a/examples/SubgridTests/StellarEvolution/run.sh
+++ b/examples/SubgridTests/StellarEvolution/run.sh
@@ -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.log
+../../swift --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.log
+../../swift --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.log
+../../swift --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.log
+../../swift --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.log
+../../swift --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
diff --git a/examples/SubgridTests/StellarEvolution/stellar_evolution.yml b/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
index 669d9e3f7c78e0ebf2b68ff0c3a3bcb5369b4a50..047a3b44027888188003921d3758d6cbab85ca7f 100644
--- a/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
+++ b/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
@@ -6,15 +6,6 @@ InternalUnitSystem:
UnitCurrent_in_cgs: 1. # Amperes
UnitTemp_in_cgs: 1. # Kelvin
-# Cosmological parameters
-Cosmology:
- h: 0.6777 # Reduced Hubble constant
- a_begin: 0.5 # Initial scale-factor of the simulation
- a_end: 1.0 # Final scale factor of the simulation
- Omega_m: 0.307 # Matter density parameter
- Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
-
# Parameters governing the time integration
TimeIntegration:
time_begin: 0 # The starting time of the simulation (in internal units).
@@ -27,6 +18,7 @@ Snapshots:
basename: stellar_evolution # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 3.e-5 # Time difference between consecutive outputs without cosmology (internal units)
+ compression: 4
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index a470eef3fafc32c99fe3a853dcf46051a3086441..f3076750478aa669bfcb557676c1eb37b72eee21 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -197,6 +197,16 @@ __attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
sp->chemistry_data.metal_mass_fraction[elem] =
data->initial_metal_mass_fraction[elem];
}
+
+ /* Initialize mass fractions for total metals and each metal individually */
+ if (data->initial_metal_mass_fraction_total != -1) {
+ sp->chemistry_data.smoothed_metal_mass_fraction_total =
+ data->initial_metal_mass_fraction_total;
+
+ for (int elem = 0; elem < chemistry_element_count; ++elem)
+ sp->chemistry_data.smoothed_metal_mass_fraction[elem] =
+ data->initial_metal_mass_fraction[elem];
+ }
}
/**
diff --git a/src/scheduler.c b/src/scheduler.c
index ac00949a9bcdb3974738215c173fdbdf7212e756..d060d5d469389b42c865653c0892472baef062f8 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -568,6 +568,7 @@ static void scheduler_splittask_hydro(struct task *t, struct scheduler *s) {
if (!is_self && !is_pair) {
t->type = task_type_none;
t->subtype = task_subtype_none;
+ t->ci = NULL;
t->cj = NULL;
t->skip = 1;
break;
@@ -767,6 +768,7 @@ static void scheduler_splittask_gravity(struct task *t, struct scheduler *s) {
if ((t->ci == NULL) || (t->type == task_type_pair && t->cj == NULL)) {
t->type = task_type_none;
t->subtype = task_subtype_none;
+ t->ci = NULL;
t->cj = NULL;
t->skip = 1;
break;
@@ -912,6 +914,7 @@ static void scheduler_splittask_fof(struct task *t, struct scheduler *s) {
t->ci->grav.count == 0 || (t->cj != NULL && t->cj->grav.count == 0)) {
t->type = task_type_none;
t->subtype = task_subtype_none;
+ t->ci = NULL;
t->cj = NULL;
t->skip = 1;
break;