Skip to content
Snippets Groups Projects
Commit f31a4aea authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Use a smaller glass file for the stellar evolution test. Also plot the stellar mass evolution.

parent 7550199d
No related branches found
No related tags found
1 merge request!787Eagle stellar evolution matthieu
#!/bin/bash #!/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
...@@ -45,7 +45,7 @@ boxsize_cgs = 1.0e1*kpc_in_cm ...@@ -45,7 +45,7 @@ boxsize_cgs = 1.0e1*kpc_in_cm
vol_cgs = boxsize_cgs**3 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 # Read particle positions and h from the glass
pos = glass["/PartType0/Coordinates"][:,:] pos = glass["/PartType0/Coordinates"][:,:]
......
...@@ -124,22 +124,26 @@ for line in eagle_data: ...@@ -124,22 +124,26 @@ for line in eagle_data:
i += 1 i += 1
# Declare arrays to store SWIFT data # Declare arrays to store SWIFT data
swift_box_mass = zeros(n_snapshots) swift_box_gas_mass = zeros(n_snapshots)
swift_box_metal_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_element_mass = zeros((n_snapshots,n_elements))
t = zeros(n_snapshots) t = zeros(n_snapshots)
# Read data from snapshots # Read data from snapshots
for i in range(n_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") sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
t[i] = sim["/Header"].attrs["Time"][0] t[i] = sim["/Header"].attrs["Time"][0]
masses = sim["/PartType0/Masses"][:] 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"][:] 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"][:][:] element_abundances = sim["/PartType0/ElementAbundance"][:][:]
for j in range(n_elements): for j in range(n_elements):
...@@ -150,25 +154,25 @@ for i in range(n_snapshots): ...@@ -150,25 +154,25 @@ for i in range(n_snapshots):
# Plot the interesting quantities # Plot the interesting quantities
figure() figure()
# Box mass -------------------------------- # Box gas mass --------------------------------
subplot(221) 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') plot(eagle_time_Gyr[1:],eagle_total_mass[:-1],linewidth=0.5,color='r',label='eagle test total')
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0) xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in total gas particle mass (Msun)", labelpad=2) ylabel("Change in total gas particle mass (Msun)", labelpad=2)
ticklabel_format(style='sci', axis='y', scilimits=(0,0)) ticklabel_format(style='sci', axis='y', scilimits=(0,0))
legend() legend()
# Box metal mass -------------------------------- # Box star mass --------------------------------
subplot(222) 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(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_metal_mass[:-1],linewidth=0.5,color='r',label='eagle test') plot(eagle_time_Gyr[1:],-eagle_total_mass[:-1],linewidth=0.5,color='r',label='eagle test total')
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0) 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)) 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'] colours = ['k','r','g','b','c','y','m','skyblue','plum']
element_names = ['H','He','C','N','O','Ne','Mg','Si','Fe'] element_names = ['H','He','C','N','O','Ne','Mg','Si','Fe']
subplot(223) subplot(223)
...@@ -177,10 +181,17 @@ for j in range(n_elements): ...@@ -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='--') plot(eagle_time_Gyr[1:],eagle_total_element_mass[:-1,j],linewidth=1,color=colours[j],linestyle='--')
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0) xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in element mass of gas particles (Msun)", labelpad=2) ylabel("Change in element mass of gas particles (Msun)", labelpad=2)
#ticklabel_format(style='sci', axis='y', scilimits=(0,0))
xscale("log") xscale("log")
yscale("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) savefig("box_evolution.png", dpi=200)
......
#!/bin/bash #!/bin/bash
# Generate the initial conditions if they are not present. # Generate the initial conditions if they are not present.
if [ ! -e glassCube_64.hdf5 ] if [ ! -e glassCube_32.hdf5 ]
then then
echo "Fetching initial glass file for the Supernovae feedback example..." echo "Fetching initial glass file for the Supernovae feedback example..."
./getGlass.sh ./getGlass.sh
...@@ -19,29 +19,8 @@ then ...@@ -19,29 +19,8 @@ then
./getYieldTable.sh ./getYieldTable.sh
fi 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 #python check_stellar_evolution.py
# python plot_box_evolution.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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment