diff --git a/examples/SubgridTests/StellarEvolution/getGlass.sh b/examples/SubgridTests/StellarEvolution/getGlass.sh
index d5c5f590ac37c9c9431d626a2ea61b0c12c1513c..ffd92e88deae6e91237059adac2a6c2067caee46 100755
--- a/examples/SubgridTests/StellarEvolution/getGlass.sh
+++ b/examples/SubgridTests/StellarEvolution/getGlass.sh
@@ -1,2 +1,2 @@
#!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
diff --git a/examples/SubgridTests/StellarEvolution/makeIC.py b/examples/SubgridTests/StellarEvolution/makeIC.py
index 469e6ef723aa6b8128e6d20bdc508a7735cbb4fb..880e47844da90d42af1f287abe7a375ad4e4639f 100644
--- a/examples/SubgridTests/StellarEvolution/makeIC.py
+++ b/examples/SubgridTests/StellarEvolution/makeIC.py
@@ -45,7 +45,7 @@ boxsize_cgs = 1.0e1*kpc_in_cm
vol_cgs = boxsize_cgs**3
#---------------------------------------------------
-glass = h5py.File("glassCube_64.hdf5", "r")
+glass = h5py.File("glassCube_32.hdf5", "r")
# Read particle positions and h from the glass
pos = glass["/PartType0/Coordinates"][:,:]
diff --git a/examples/SubgridTests/StellarEvolution/plot_box_evolution.py b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
index d72c50c1f2152d4201dacb178f6bfe9e5a624831..15fa363ca8e606f1d842b5703f62955b141df987 100644
--- a/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
+++ b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
@@ -124,22 +124,26 @@ for line in eagle_data:
i += 1
# Declare arrays to store SWIFT data
-swift_box_mass = zeros(n_snapshots)
-swift_box_metal_mass = zeros(n_snapshots)
+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))
t = zeros(n_snapshots)
# Read data from snapshots
for i in range(n_snapshots):
- print("reading snapshot "+str(i))
+ #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_mass[i] = np.sum(masses)
-
+ swift_box_gas_mass[i] = np.sum(masses)
+
+ star_masses = sim["/PartType4/Masses"][:]
+ swift_box_star_mass[i] = np.sum(star_masses)
+
metallicities = sim["/PartType0/Metallicity"][:]
- swift_box_metal_mass[i] = np.sum(metallicities * masses)
+ swift_box_gas_metal_mass[i] = np.sum(metallicities * masses)
element_abundances = sim["/PartType0/ElementAbundance"][:][:]
for j in range(n_elements):
@@ -150,25 +154,25 @@ for i in range(n_snapshots):
# Plot the interesting quantities
figure()
-# Box mass --------------------------------
+# Box gas mass --------------------------------
subplot(221)
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (swift_box_mass[1:] - swift_box_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', marker = "*", ms=0.5, label='swift')
+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')
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 metal mass --------------------------------
+# Box star mass --------------------------------
subplot(222)
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (swift_box_metal_mass[1:] - swift_box_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')
+plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (swift_box_star_mass[1:] - swift_box_star_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')
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
-ylabel("Change in total metal mass of gas particles (Msun)", labelpad=2)
+ylabel("Change in total star particle mass (Msun)", labelpad=2)
ticklabel_format(style='sci', axis='y', scilimits=(0,0))
-#legend()
+legend()
-# Box element mass --------------------------------
+# 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)
@@ -177,10 +181,17 @@ for j in range(n_elements):
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)
-#ticklabel_format(style='sci', axis='y', scilimits=(0,0))
xscale("log")
yscale("log")
-legend(bbox_to_anchor=(1.1, 1.))
+legend(bbox_to_anchor=(2.1, 0.5), ncol=3)
+
+# 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.png", dpi=200)
diff --git a/examples/SubgridTests/StellarEvolution/run.sh b/examples/SubgridTests/StellarEvolution/run.sh
index 2aa980b582be1068a7c3f80d06062855e346e46d..8fa1e0b56ec39552c1c9f60c40290c551276b7c1 100755
--- a/examples/SubgridTests/StellarEvolution/run.sh
+++ b/examples/SubgridTests/StellarEvolution/run.sh
@@ -1,7 +1,7 @@
#!/bin/bash
# Generate the initial conditions if they are not present.
-if [ ! -e glassCube_64.hdf5 ]
+if [ ! -e glassCube_32.hdf5 ]
then
echo "Fetching initial glass file for the Supernovae feedback example..."
./getGlass.sh
@@ -19,29 +19,8 @@ then
./getYieldTable.sh
fi
-# Run mass enrichment check
-#../../swift --limiter --feedback --stars --hydro --external-gravity --threads=4 --param=Stars:energy_testing:0 --param=TimeIntegration:time_end:1e-4 --param=Stars:feedback_timescale:1e-4 stellar_evolution.yml 2>&1 | tee output_enrichment.log
-#
-#python check_stellar_evolution.py
-#
-## Run continuous heating check
-#../../swift --limiter --feedback --stars --hydro --external-gravity --threads=4 --param=Stars:energy_testing:1 --param=Stars:continuous_heating:1 stellar_evolution.yml 2>&1 | tee output_continuous.log
-#
-#python check_continuous_heating.py
-#
-## Run stochastic check
-#../../swift --limiter --feedback --stars --hydro --external-gravity --threads=4 --param=Stars:energy_testing:1 stellar_evolution.yml 2>&1 | tee output_stochastic_1.log
-#
-#python check_stochastic_heating.py
-../../swift --limiter --feedback --stars --hydro --external-gravity --threads=4 --param=Stars:energy_testing:1 stellar_evolution.yml 2>&1 | tee output.log
+../../swift --limiter --feedback --stars --hydro --external-gravity --threads=4 stellar_evolution.yml 2>&1 | tee output.log
-python check_stochastic_heating.py
-#
-#../../swift --limiter --feedback --stars --hydro --external-gravity --threads=4 --param=Stars:energy_testing:1 stellar_evolution.yml 2>&1 | tee output_stochastic_3.log
-#
-#python check_stochastic_heating.py
-#
-#../../swift --limiter --feedback --stars --hydro --external-gravity --threads=4 --param=Stars:energy_testing:1 stellar_evolution.yml 2>&1 | tee output_stochastic_4.log
-#
-#python check_stochastic_heating.py
+#python check_stellar_evolution.py
+python plot_box_evolution.py