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/configure.ac b/configure.ac
index fca63462f69779b03663ad30b23ce876c6c85002..885b93479cf333c7cac03b74e18983c6aec6a07f 100644
--- a/configure.ac
+++ b/configure.ac
@@ -432,14 +432,8 @@ HAVEVECTORIZATION=0
if test "$enable_opt" = "yes" ; then
- # Add code optimisation flags and tuning to host. This is a funny macro
- # that does not like CFLAGS being already set. Work around that as we have
- # at least set it to "", so it is set.
- ac_test_CFLAGS="no"
- old_CFLAGS="$CFLAGS"
+ # Choose the best flags for this compiler and architecture
AX_CC_MAXOPT
- ac_test_CFLAGS="yes"
- CFLAGS="$old_CFLAGS $CFLAGS"
# Check SSE & AVX support (some overlap with AX_CC_MAXOPT).
# Don't use the SIMD_FLAGS result with Intel compilers. The -x<code>
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/plot_box_evolution.py b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
index bcfa85a1afac021e280a797641dd677bccd291d3..96f4acf4352660d9ff05ffa62b83b156e654da9f 100644
--- a/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
+++ b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
@@ -104,15 +104,15 @@ 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/Metallicity"][0]
+ Z_star = sim["/PartType4/MetalMassFractions"][0]
star_masses = sim["/PartType4/Masses"][:]
swift_box_star_mass[i] = np.sum(star_masses)
- metallicities = sim["/PartType0/Metallicities"][:]
+ metallicities = sim["/PartType0/MetalMassFractions"][:]
swift_box_gas_metal_mass[i] = np.sum(metallicities * masses)
element_abundances = sim["/PartType0/ElementMassFractions"][:][:]
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/examples/nIFTyCluster/Baryonic/convert_to_swift.py b/examples/nIFTyCluster/Baryonic/convert_to_swift.py
index 1de42ac8021e5fa05c4bfe78df8a303669db4a6e..323ca05f1c1df5bbe2744fb882be38d0834c82e3 100644
--- a/examples/nIFTyCluster/Baryonic/convert_to_swift.py
+++ b/examples/nIFTyCluster/Baryonic/convert_to_swift.py
@@ -18,12 +18,12 @@ filename = "/cosma7/data/dp004/jlvc76/nIFTy/IC_CLUSTER_00019"
length = unyt.kpc
mass = 1e10 * unyt.msun
-time = unyt.s * unyt.kpc / unyt.km
+time = (1.0 * unyt.s * unyt.kpc / unyt.km).to("s")
velocity = length / time
energy_per_unit_mass = (length / time) ** 2
-nifty_units = unyt.UnitSystem("nifty", length, mass, time)
+nifty_units = unyt.UnitSystem("nifty", 1e3 * length, mass, time)
writer = Writer(
unit_system=nifty_units,
diff --git a/examples/nIFTyCluster/Baryonic/nifty.yml b/examples/nIFTyCluster/Baryonic/nifty.yml
index dcf1f0c1e399048984d9c397f4638e75b8f86cf2..d520d5bfba6ca466c30a0bd234a5457d7edd8fdf 100644
--- a/examples/nIFTyCluster/Baryonic/nifty.yml
+++ b/examples/nIFTyCluster/Baryonic/nifty.yml
@@ -1,7 +1,7 @@
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
- UnitLength_in_cgs: 3.08567758e21 # kpc in centimeters
+ 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
@@ -46,10 +46,9 @@ Gravity:
softening_ratio_background: 0.04 # 1/25th of mean inter-particle separation
mesh_side_length: 512
-
# Parameters for the hydrodynamics scheme
SPH:
- resolution_eta: 1.48691 # 100 ngb with wendland0-c2
+ resolution_eta: 1.2 # ~50 ngb with Wendland-C2
h_min_ratio: 0.1 # Minimal smoothing in units of softening.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
minimal_temperature: 100 # (internal units)
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index bee8608820f092bf4534a192c7d3dab65a603702..7d8a5886ebcde64a39defa24b9f1e9fbd2232642 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -48,6 +48,13 @@ SPH:
diffusion_alpha_max: 1.0 # (Optional) Override the maximal thermal diffusion coefficient that is allowed for a given particle.
diffusion_alpha_min: 0.0 # (Optional) Override the minimal thermal diffusion coefficient that is allowed for a given particle.
+# Parameters of the stars neighbour search
+Stars:
+ resolution_eta: 1.2348 # (Optional) Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). Defaults to the SPH value.
+ h_tolerance: 1e-4 # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths. Defaults to the SPH value.
+ max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length. Defaults to the SPH value.
+ max_volume_change: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step. Defaults to the SPH value.
+ birth_time: -1 # (Optional) Initial birth time of *all* the stars. If not -1, this value will overwrite all the values read from the ICs.
# Parameters for the self-gravity scheme
Gravity:
diff --git a/m4/ax_cc_maxopt.m4 b/m4/ax_cc_maxopt.m4
index 7523d9d09ee741914fe9fe1cb22d3e0550763233..72407f3b6a5f29758fbb7b41afda9357ce4210ec 100644
--- a/m4/ax_cc_maxopt.m4
+++ b/m4/ax_cc_maxopt.m4
@@ -1,5 +1,5 @@
# ===========================================================================
-# http://www.gnu.org/software/autoconf-archive/ax_cc_maxopt.html
+# https://www.gnu.org/software/autoconf-archive/ax_cc_maxopt.html
# ===========================================================================
#
# SYNOPSIS
@@ -40,7 +40,7 @@
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
-# with this program. If not, see <http://www.gnu.org/licenses/>.
+# with this program. If not, see <https://www.gnu.org/licenses/>.
#
# As a special exception, the respective Autoconf Macro's copyright owner
# gives unlimited permission to copy, distribute and modify the configure
@@ -55,7 +55,7 @@
# modified version of the Autoconf Macro, you may extend this special
# exception to the GPL to apply to your modified version as well.
-#serial 16
+#serial 18
AC_DEFUN([AX_CC_MAXOPT],
[
@@ -68,19 +68,18 @@ AC_ARG_ENABLE(portable-binary, [AS_HELP_STRING([--enable-portable-binary], [disa
# Try to determine "good" native compiler flags if none specified via CFLAGS
if test "$ac_test_CFLAGS" != "set"; then
- CFLAGS=""
case $ax_cv_c_compiler_vendor in
- dec) CFLAGS="-newc -w0 -O5 -ansi_alias -ansi_args -fp_reorder -tune host"
+ dec) CFLAGS="$CFLAGS -newc -w0 -O5 -ansi_alias -ansi_args -fp_reorder -tune host"
if test "x$acx_maxopt_portable" = xno; then
CFLAGS="$CFLAGS -arch host"
fi;;
- sun) CFLAGS="-native -fast -xO5 -dalign"
+ sun) CFLAGS="$CFLAGS -native -fast -xO5 -dalign"
if test "x$acx_maxopt_portable" = xyes; then
CFLAGS="$CFLAGS -xarch=generic"
fi;;
- hp) CFLAGS="+Oall +Optrs_ansi +DSnative"
+ hp) CFLAGS="$CFLAGS +Oall +Optrs_ansi +DSnative"
if test "x$acx_maxopt_portable" = xyes; then
CFLAGS="$CFLAGS +DAportable"
fi;;
@@ -91,8 +90,8 @@ if test "$ac_test_CFLAGS" != "set"; then
xlc_opt="-qtune=auto"
fi
AX_CHECK_COMPILE_FLAG($xlc_opt,
- CFLAGS="-O3 -qansialias -w $xlc_opt",
- [CFLAGS="-O3 -qansialias -w"
+ CFLAGS="$CFLAGS -O3 -qansialias -w $xlc_opt",
+ [CFLAGS="$CFLAGS -O3 -qansialias -w"
echo "******************************************************"
echo "* You seem to have the IBM C compiler. It is *"
echo "* recommended for best performance that you use: *"
@@ -105,7 +104,7 @@ if test "$ac_test_CFLAGS" != "set"; then
echo "******************************************************"])
;;
- intel) CFLAGS="-O3 -ansi-alias"
+ intel) CFLAGS="$CFLAGS -O3 -ansi-alias"
if test "x$acx_maxopt_portable" = xno; then
icc_archflag=unknown
icc_flags=""
@@ -151,7 +150,7 @@ if test "$ac_test_CFLAGS" != "set"; then
clang)
# default optimization flags for clang on all systems
- CFLAGS="-O3 -fomit-frame-pointer"
+ CFLAGS="$CFLAGS -O3 -fomit-frame-pointer"
# Always good optimisation to have
AX_CHECK_COMPILE_FLAG(-fstrict-aliasing, CFLAGS="$CFLAGS -fstrict-aliasing")
@@ -167,7 +166,7 @@ if test "$ac_test_CFLAGS" != "set"; then
gnu)
# default optimization flags for gcc on all systems
- CFLAGS="-O3 -fomit-frame-pointer"
+ CFLAGS="$CFLAGS -O3 -fomit-frame-pointer"
# -malign-double for x86 systems
AX_CHECK_COMPILE_FLAG(-malign-double, CFLAGS="$CFLAGS -malign-double")
@@ -187,7 +186,7 @@ if test "$ac_test_CFLAGS" != "set"; then
microsoft)
# default optimization flags for MSVC opt builds
- CFLAGS="-O2"
+ CFLAGS="$CFLAGS -O2"
;;
esac
@@ -199,7 +198,7 @@ if test "$ac_test_CFLAGS" != "set"; then
echo "* (otherwise, a default of CFLAGS=-O3 will be used) *"
echo "********************************************************"
echo ""
- CFLAGS="-O3"
+ CFLAGS="$CFLAGS -O3"
fi
AX_CHECK_COMPILE_FLAG($CFLAGS, [], [
@@ -210,7 +209,6 @@ if test "$ac_test_CFLAGS" != "set"; then
echo "* Use ./configure CFLAGS=... to specify your own flags *"
echo "********************************************************"
echo ""
- CFLAGS=""
])
fi
diff --git a/src/cell.c b/src/cell.c
index 5191cde2facc209df16a8e42788655611cebcb15..4b9746e92e31f2a5ae6c4b5d01ade8d83939e412 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -3508,7 +3508,9 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
if (c->timestep != NULL) scheduler_activate(s, c->timestep);
if (c->hydro.end_force != NULL) scheduler_activate(s, c->hydro.end_force);
if (c->hydro.cooling != NULL) scheduler_activate(s, c->hydro.cooling);
+#ifdef WITH_LOGGER
if (c->logger != NULL) scheduler_activate(s, c->logger);
+#endif
if (c->top->hydro.star_formation != NULL) {
cell_activate_star_formation_tasks(c->top, s);
@@ -3663,7 +3665,9 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
if (c->grav.mesh != NULL) scheduler_activate(s, c->grav.mesh);
if (c->grav.long_range != NULL) scheduler_activate(s, c->grav.long_range);
if (c->grav.end_force != NULL) scheduler_activate(s, c->grav.end_force);
+#ifdef WITH_LOGGER
if (c->logger != NULL) scheduler_activate(s, c->logger);
+#endif
}
return rebuild;
@@ -3905,7 +3909,9 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s,
if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
if (c->timestep != NULL) scheduler_activate(s, c->timestep);
+#ifdef WITH_LOGGER
if (c->logger != NULL) scheduler_activate(s, c->logger);
+#endif
}
}
@@ -4183,7 +4189,9 @@ int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s) {
if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
if (c->timestep != NULL) scheduler_activate(s, c->timestep);
+#ifdef WITH_LOGGER
if (c->logger != NULL) scheduler_activate(s, c->logger);
+#endif
}
return rebuild;
diff --git a/src/cell.h b/src/cell.h
index 3bdae3395452df21455068fa8e7a2707749f067b..8067a3189818ab8738de848ea698fbf25c78ebba 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -755,8 +755,10 @@ struct cell {
/*! The task to limit the time-step of inactive particles */
struct task *timestep_limiter;
+#ifdef WITH_LOGGER
/*! The logger task */
struct task *logger;
+#endif
/*! Minimum dimension, i.e. smallest edge of this cell (min(width)). */
float dmin;
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/engine_maketasks.c b/src/engine_maketasks.c
index 6a8341068cb87d2994a8d658a9a7255462ef294e..05bde9091dd55904063133c1a70cc004f0a05512 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -789,24 +789,33 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
c->kick1 = scheduler_addtask(s, task_type_kick1, task_subtype_none, 0, 0,
c, NULL);
+ c->kick2 = scheduler_addtask(s, task_type_kick2, task_subtype_none, 0, 0,
+ c, NULL);
+
#if defined(WITH_LOGGER)
+ /* Add the hydro logger task. */
c->logger = scheduler_addtask(s, task_type_logger, task_subtype_none, 0,
0, c, NULL);
-#endif
- c->kick2 = scheduler_addtask(s, task_type_kick2, task_subtype_none, 0, 0,
- c, NULL);
+ /* Add the kick2 dependency */
+ scheduler_addunlock(s, c->kick2, c->logger);
+
+ /* Create a variable in order to avoid to many ifdef */
+ struct task *kick2_or_logger = c->logger;
+#else
+ struct task *kick2_or_logger = c->kick2;
+#endif
/* Add the time-step calculation task and its dependency */
c->timestep = scheduler_addtask(s, task_type_timestep, task_subtype_none,
0, 0, c, NULL);
- scheduler_addunlock(s, c->kick2, c->timestep);
+ scheduler_addunlock(s, kick2_or_logger, c->timestep);
scheduler_addunlock(s, c->timestep, c->kick1);
/* Subgrid tasks: star formation */
if (with_star_formation && c->hydro.count > 0) {
- scheduler_addunlock(s, c->kick2, c->top->hydro.star_formation);
+ scheduler_addunlock(s, kick2_or_logger, c->top->hydro.star_formation);
scheduler_addunlock(s, c->top->hydro.star_formation, c->timestep);
}
@@ -819,10 +828,6 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
scheduler_addunlock(s, c->timestep, c->timestep_limiter);
scheduler_addunlock(s, c->timestep_limiter, c->kick1);
}
-
-#if defined(WITH_LOGGER)
- scheduler_addunlock(s, c->kick1, c->logger);
-#endif
}
} else { /* We are above the super-cell so need to go deeper */
@@ -1094,7 +1099,11 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
c->stars.ghost = scheduler_addtask(s, task_type_stars_ghost,
task_subtype_none, 0, 0, c, NULL);
+#ifdef WITH_LOGGER
+ scheduler_addunlock(s, c->super->logger, c->stars.stars_in);
+#else
scheduler_addunlock(s, c->super->kick2, c->stars.stars_in);
+#endif
scheduler_addunlock(s, c->stars.stars_out, c->super->timestep);
if (with_star_formation && c->hydro.count > 0) {
@@ -1124,7 +1133,11 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
c->black_holes.swallow_ghost[2] = scheduler_addtask(
s, task_type_bh_swallow_ghost3, task_subtype_none, 0, 0, c, NULL);
+#ifdef WITH_LOGGER
+ scheduler_addunlock(s, c->super->logger, c->black_holes.black_holes_in);
+#else
scheduler_addunlock(s, c->super->kick2, c->black_holes.black_holes_in);
+#endif
scheduler_addunlock(s, c->black_holes.black_holes_out,
c->super->timestep);
}
@@ -1758,8 +1771,21 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
const enum task_types t_type = t->type;
const enum task_subtypes t_subtype = t->subtype;
const long long flags = t->flags;
- struct cell *ci = t->ci;
- struct cell *cj = t->cj;
+ struct cell *const ci = t->ci;
+ struct cell *const cj = t->cj;
+
+ /* Escape early */
+ if (t->type == task_type_none) continue;
+
+#ifdef WITH_LOGGER
+ struct task *const ci_super_kick2_or_logger = ci->super->logger;
+ struct task *const cj_super_kick2_or_logger =
+ (cj == NULL) ? NULL : cj->super->logger;
+#else
+ struct task *const ci_super_kick2_or_logger = ci->super->kick2;
+ struct task *const cj_super_kick2_or_logger =
+ (cj == NULL) ? NULL : cj->super->kick2;
+#endif
/* Sort tasks depend on the drift of the cell (gas version). */
if (t_type == task_type_sort && ci->nodeID == nodeID) {
@@ -1907,7 +1933,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_limiter) {
- scheduler_addunlock(sched, ci->super->kick2, t_limiter);
+ scheduler_addunlock(sched, ci_super_kick2_or_logger, t_limiter);
scheduler_addunlock(sched, t_limiter, ci->super->timestep);
scheduler_addunlock(sched, t_limiter, ci->super->timestep_limiter);
}
@@ -2097,7 +2123,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_limiter) {
- scheduler_addunlock(sched, ci->super->kick2, t_limiter);
+ scheduler_addunlock(sched, ci_super_kick2_or_logger, t_limiter);
scheduler_addunlock(sched, t_limiter, ci->super->timestep);
scheduler_addunlock(sched, t_limiter, ci->super->timestep_limiter);
}
@@ -2175,7 +2201,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_limiter) {
- scheduler_addunlock(sched, cj->super->kick2, t_limiter);
+ scheduler_addunlock(sched, cj_super_kick2_or_logger, t_limiter);
scheduler_addunlock(sched, t_limiter, cj->super->timestep);
scheduler_addunlock(sched, t_limiter, cj->super->timestep_limiter);
}
@@ -2343,7 +2369,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_limiter) {
- scheduler_addunlock(sched, ci->super->kick2, t_limiter);
+ scheduler_addunlock(sched, ci_super_kick2_or_logger, t_limiter);
scheduler_addunlock(sched, t_limiter, ci->super->timestep);
scheduler_addunlock(sched, t_limiter, ci->super->timestep_limiter);
}
@@ -2537,7 +2563,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_limiter) {
- scheduler_addunlock(sched, ci->super->kick2, t_limiter);
+ scheduler_addunlock(sched, ci_super_kick2_or_logger, t_limiter);
scheduler_addunlock(sched, t_limiter, ci->super->timestep);
scheduler_addunlock(sched, t_limiter, ci->super->timestep_limiter);
}
@@ -2617,7 +2643,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_limiter) {
- scheduler_addunlock(sched, cj->super->kick2, t_limiter);
+ scheduler_addunlock(sched, cj_super_kick2_or_logger, t_limiter);
scheduler_addunlock(sched, t_limiter, cj->super->timestep);
scheduler_addunlock(sched, t_limiter, cj->super->timestep_limiter);
}
@@ -2962,7 +2988,7 @@ void engine_make_fof_tasks(struct engine *e) {
tic = getticks();
/* Split the tasks. */
- scheduler_splittasks(sched, /*fof_tasks=*/1);
+ scheduler_splittasks(sched, /*fof_tasks=*/1, e->verbose);
if (e->verbose)
message("Splitting FOF tasks took %.3f %s.",
@@ -3038,7 +3064,7 @@ void engine_maketasks(struct engine *e) {
tic2 = getticks();
/* Split the tasks. */
- scheduler_splittasks(sched, /*fof_tasks=*/0);
+ scheduler_splittasks(sched, /*fof_tasks=*/0, e->verbose);
if (e->verbose)
message("Splitting tasks took %.3f %s.",
diff --git a/src/runner.c b/src/runner.c
index e1c3b5139d6cad9095f46f10b47584daffc86892..db7e512873b51a7329e19e75a763b69521efb0eb 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -4934,7 +4934,7 @@ void runner_do_logger(struct runner *r, struct cell *c, int timer) {
const int count = c->hydro.count;
/* Anything to do here? */
- if (!cell_is_starting_hydro(c, e) && !cell_is_starting_gravity(c, e)) return;
+ if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e)) return;
/* Recurse? Avoid spending too much time in useless cells. */
if (c->split) {
@@ -4952,7 +4952,7 @@ void runner_do_logger(struct runner *r, struct cell *c, int timer) {
/* If particle needs to be log */
/* This is the same function than part_is_active, except for
* debugging checks */
- if (part_is_starting(p, e)) {
+ if (part_is_active(p, e)) {
if (logger_should_write(&xp->logger_data, e->logger)) {
/* Write particle */
diff --git a/src/scheduler.c b/src/scheduler.c
index ac00949a9bcdb3974738215c173fdbdf7212e756..a1c883fabdac3397692ac7bb93dbd5586fb3f19b 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;
@@ -1026,7 +1029,17 @@ void scheduler_splittasks_mapper(void *map_data, int num_elements,
* @param fof_tasks Are we splitting the FOF tasks (1)? Or the regular tasks
* (0)?
*/
-void scheduler_splittasks(struct scheduler *s, const int fof_tasks) {
+void scheduler_splittasks(struct scheduler *s, const int fof_tasks,
+ const int verbose) {
+
+ if (verbose) {
+ message("space_subsize_self_hydro= %d", space_subsize_self_hydro);
+ message("space_subsize_pair_hydro= %d", space_subsize_pair_hydro);
+ message("space_subsize_self_stars= %d", space_subsize_self_stars);
+ message("space_subsize_pair_stars= %d", space_subsize_pair_stars);
+ message("space_subsize_self_grav= %d", space_subsize_self_grav);
+ message("space_subsize_pair_grav= %d", space_subsize_pair_grav);
+ }
if (fof_tasks) {
/* Call the mapper on each current task. */
diff --git a/src/scheduler.h b/src/scheduler.h
index ac9bc754db1ea03f91c0ce522dc325c18d25ec09..694a88017b863d29ad9d15a4ae7d5acac5cf3223 100644
--- a/src/scheduler.h
+++ b/src/scheduler.h
@@ -190,7 +190,8 @@ void scheduler_reweight(struct scheduler *s, int verbose);
struct task *scheduler_addtask(struct scheduler *s, enum task_types type,
enum task_subtypes subtype, int flags,
int implicit, struct cell *ci, struct cell *cj);
-void scheduler_splittasks(struct scheduler *s, const int fof_tasks);
+void scheduler_splittasks(struct scheduler *s, const int fof_tasks,
+ const int verbose);
struct task *scheduler_done(struct scheduler *s, struct task *t);
struct task *scheduler_unlock(struct scheduler *s, struct task *t);
void scheduler_addunlock(struct scheduler *s, struct task *ta, struct task *tb);
diff --git a/src/stars/EAGLE/stars.h b/src/stars/EAGLE/stars.h
index 09926d6a8bdd5af4eff120a9a6c57d15eb7e49c4..f102ddb22c0075d33b02c726d0447b64fa9d3df7 100644
--- a/src/stars/EAGLE/stars.h
+++ b/src/stars/EAGLE/stars.h
@@ -65,7 +65,8 @@ __attribute__((always_inline)) INLINE static void stars_first_init_spart(
sp->time_bin = 0;
sp->birth_density = 0.f;
sp->f_E = -1.f;
- sp->birth_time = stars_properties->spart_first_init_birth_time;
+ if (stars_properties->spart_first_init_birth_time != -1.f)
+ sp->birth_time = stars_properties->spart_first_init_birth_time;
stars_init_spart(sp);
}
diff --git a/src/stars/EAGLE/stars_io.h b/src/stars/EAGLE/stars_io.h
index b67c57078796a4be3d34d828884cedcff808614e..b91b5cf94595a05acd280cfd4f51755f91cce04d 100644
--- a/src/stars/EAGLE/stars_io.h
+++ b/src/stars/EAGLE/stars_io.h
@@ -35,7 +35,7 @@ INLINE static void stars_read_particles(struct spart *sparts,
int *num_fields) {
/* Say how much we want to read */
- *num_fields = 6;
+ *num_fields = 7;
/* List what we want to read */
list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
@@ -50,6 +50,8 @@ INLINE static void stars_read_particles(struct spart *sparts,
UNIT_CONV_LENGTH, sparts, h);
list[5] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
sparts, mass_init);
+ list[6] = io_make_input_field("StellarFormationTime", FLOAT, 1, OPTIONAL,
+ UNIT_CONV_NO_UNITS, sparts, birth_time);
}
INLINE static void convert_spart_pos(const struct engine *e,
@@ -218,7 +220,7 @@ INLINE static void stars_props_init(struct stars_props *sp,
/* Read birth time to set all stars in ICs to (defaults to -1 to indicate star
* present in ICs) */
sp->spart_first_init_birth_time =
- parser_get_opt_param_float(params, "Stars:birth_time", -1);
+ parser_get_opt_param_float(params, "Stars:birth_time", -1.f);
}
/**
diff --git a/src/stars/EAGLE/stars_part.h b/src/stars/EAGLE/stars_part.h
index bed54bb175026ae930668fb102dcfef9c9af76d7..4502b10edb7b646e9ba845e1ffffbb9255cdc01c 100644
--- a/src/stars/EAGLE/stars_part.h
+++ b/src/stars/EAGLE/stars_part.h
@@ -153,7 +153,7 @@ struct stars_props {
/*! Maximal change of h over one time-step */
float log_max_h_change;
- /*! Value to set birth time of stars read from ICs */
+ /*! Value to set birth time of stars read from ICs if not set to -1 */
float spart_first_init_birth_time;
};