diff --git a/examples/SubgridTests/StellarEvolution/check_continuous_heating.py b/examples/SubgridTests/StellarEvolution/check_continuous_heating.py
index 5035763c1b8167384a4302cefb37a13169d07f62..1a998d3d9f25d5620755d160c33d3a49342ce1fc 100644
--- a/examples/SubgridTests/StellarEvolution/check_continuous_heating.py
+++ b/examples/SubgridTests/StellarEvolution/check_continuous_heating.py
@@ -90,7 +90,6 @@ smoothing_length_sparts = zeros(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))
diff --git a/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py b/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py
index c0426f31824bd7d27b3f8f28833e55cdbf350556..da837540041a9295a33b55e16b5e996394576cd7 100644
--- a/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py
+++ b/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py
@@ -1,3 +1,10 @@
+# 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 *
@@ -6,15 +13,6 @@ import os.path
import numpy as np
import glob
-
-# Function to determine distance between part and spart
-# p: part coordinates
-# s: spart coordinates
-def distance(p,s):
- dist2 = (p[0] - s[0]) * (p[0] - s[0]) + (p[1] - s[1]) * (p[1] - s[1]) +(p[2] - s[2]) * (p[2] - s[2])
- return sqrt(dist2)
-
-
# 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
@@ -72,6 +70,12 @@ 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]
@@ -79,105 +83,55 @@ 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,n_snapshots))
velocity_parts = zeros((n_parts,3,n_snapshots))
-speed_parts = zeros((n_parts,n_snapshots))
-coord_sparts = zeros((3,n_snapshots))
-smoothing_length_parts = zeros((n_parts,n_snapshots))
-distances = zeros((n_parts,n_snapshots))
-smoothing_length_sparts = zeros(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):
+#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))
- abundances[:,:,i] = sim["/PartType0/ElementAbundance"]
- metallicity[:,i] = sim["/PartType0/Metallicity"]
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/InternalEnergy"]
velocity_parts[:,:,i] = sim["/PartType0/Velocities"]
- coord_parts[:,:,i] = sim["/PartType0/Coordinates"]
- coord_sparts[:,i] = sim["/PartType4/Coordinates"]
- smoothing_length_parts[:,i] = sim["/PartType0/SmoothingLength"]
- smoothing_length_sparts[i] = sim["/PartType4/SmoothingLength"][0]
time[i] = sim["/Header"].attrs["Time"][0]
# Check that the total amount of enrichment is as expected.
-# Define tolerance
-eps = 0.01
-
-#print smoothing length maximums
-#for i in range(n_snapshots):
-# max_smoothing_length_parts = np.max(smoothing_length_parts[:,i]*unit_length_in_cgs)
-# max_smoothing_length_sparts = np.max(smoothing_length_sparts[i]*unit_length_in_cgs)
-# for j in range(n_parts):
-# distances[j,i] = distance(coord_parts[j,:,i],coord_sparts[:,i])
-# min_distance_to_spart = np.min(distances[:,i])
-# print("snapshot "+ str(i) + " max smoothing length parts cgs " + str(max_smoothing_length_parts) + " max smoothing length sparts cgs " + str(max_smoothing_length_sparts) + " boxsize " + str(boxSize * unit_length_in_cgs) + " min distance to spart " + str(min_distance_to_spart))
-
+# 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 = np.sum(np.multiply(vel2,masses)*0.5,axis = 0)
-total_energy = np.sum(np.multiply(internal_energy,masses),axis = 0)
-total_energy_released = total_energy[n_snapshots-1] - total_energy[0] + total_kinetic_energy[n_snapshots-1] - total_kinetic_energy[0]
-
-# Find out how many total sn should go off in simulation time
-feedback_data = "feedback_properties.dat"
-heating_probability = 0
-energy_released = 0
-num_heating_events = 0
-with open(feedback_data) as f:
- num_sn = float(f.readline().strip())
- total_time = float(f.readline().strip())
- # Find out how many heating events occurred
- while True:
- num = f.readline().strip()
- if not num:
- break
- heating_probability += float(num.split()[0])
- energy_released += float(num.split()[1])
- num_heating_events += 1
-total_sn = num_sn * time[n_snapshots-1]/total_time
-print("total sn " + str(total_sn) + " fraction time elapsed " + str(time[n_snapshots-1]/total_time))
+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 = total_sn * energy_per_sn
-print("heating probability " + str(heating_probability) + " energy released " + str(energy_released) + " num_heating_events*energy_per_sn " + str(num_heating_events*energy_per_sn) + " expected energy released " + str(expected_energy_released))
+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 - expected_energy_released)/expected_energy_released < eps:
- print("total stochastic energy release consistent with expectation. total stochastic energy release "+str(total_energy_released)+" expected "+ str(expected_energy_released) + " initial total internal energy "+ str(total_energy[0] + total_kinetic_energy[0]))
+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)+" expected "+ str(expected_energy_released) + " initial total internal energy "+ str(total_energy[0] + total_kinetic_energy[0]) + " energy change fraction of total " + str(total_energy_released/(total_energy[0]+total_kinetic_energy[0])))
+ 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(total_energy + total_kinetic_energy,color='k')
-xlabel("snapshot")
-ylabel("total energy")
-savefig("stellar_evolution_total_energy.png", dpi=200)
+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/StellarEvolution/makeIC.py b/examples/SubgridTests/StellarEvolution/makeIC.py
index 880e47844da90d42af1f287abe7a375ad4e4639f..469e6ef723aa6b8128e6d20bdc508a7735cbb4fb 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_32.hdf5", "r")
+glass = h5py.File("glassCube_64.hdf5", "r")
# Read particle positions and h from the glass
pos = glass["/PartType0/Coordinates"][:,:]
diff --git a/examples/SubgridTests/StellarEvolution/stellar_evolution.yml b/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
index 0d73d93bd0dec69466e9d03517ea6c3d0ceca343..5eeb38055a00864e52e9bb44728c6989f1e4b1d5 100644
--- a/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
+++ b/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
@@ -20,13 +20,13 @@ TimeIntegration:
time_begin: 0 # The starting time of the simulation (in internal units).
time_end: 5e-3 # The end time of the simulation (in internal units).
dt_min: 1e-12 # The minimal time-step size of the simulation (in internal units).
- dt_max: 5e-6 # The maximal time-step size of the simulation (in internal units).
+ dt_max: 5e-7 # 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
time_first: 0. # Time of the first output (in internal units)
- delta_time: 1.e-4 # Time difference between consecutive outputs without cosmology (internal units)
+ delta_time: 1.e-5 # Time difference between consecutive outputs without cosmology (internal units)
#delta_time: 1.02 # Time difference between consecutive outputs with cosmology (internal units)
scale_factor_first: 0.5
compression: 4
@@ -34,7 +34,7 @@ Snapshots:
# Parameters governing the conserved quantities statistics
Statistics:
time_first: 0.
- delta_time: 1.e-4 # non cosmology time between statistics output
+ delta_time: 1.e-5 # non cosmology time between statistics output
#delta_time: 1.02 # cosmology time between statistics output
scale_factor_first: 0.5
@@ -91,10 +91,10 @@ EAGLEChemistry: # Solar abundances
init_abundance_Iron: 0.0
EAGLEFeedback:
- filename: ./yieldtables/
- imf_model: Chabrier
- continuous_heating_switch: 0
- SNIa_timescale_Gyr: 2.0
- SNIa_efficiency: 2.e-3
- SNII_wind_delay_Gyr: 0.03
- SNe_heating_temperature_K: 3.16228e7
+ filename: ./yieldtables/
+ imf_model: Chabrier
+ continuous_heating_switch: 0
+ SNIa_timescale_Gyr: 2.0
+ SNIa_efficiency: 2.e-3
+ SNII_wind_delay_Gyr: 0.03
+ SNe_heating_temperature_K: 3.16228e7