tidied up stochastic energy injection test

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))

examples/SubgridTests/StellarEvolution/check_stochastic_heating.py

+# 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)
 

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"][:,:]

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