In [19]:
%matplotlib inline
import numpy as np
import yt
from galaxy_analysis import Galaxy
import matplotlib.pyplot as plt
from galaxy_analysis.plot.plot_styles import *
import deepdish as dd
import glob
from galaxy_analysis.utilities import utilities
In [4]:
wdir = "/media/aemerick/Seagate Expansion Drive/FlatironData/run11_30km/sn_H2atten_H2sh/"
ds = yt.load(wdir + 'DD0239/DD0239')
data = ds.all_data()
#gal = Galaxy('DD0369',wdir=wdir)
In [5]:
#gal.ds.derived_field_list
In [6]:
data['thermal_energy']
Out[6]:
In [7]:
#
# Define region
#
dL = 400 * yt.units.pc
region = None
above_region = data.cut_region(["(obj['cylindrical_z'].in_units('kpc') > 0.8) & (obj['cylindrical_z'].in_units('kpc') < 1.2) & (obj['cylindrical_radius'].in_units('kpc') < 2)"])
above_hot_region = above_region.cut_region(["obj['temperature'].in_units('K') > 3.0E5"])
above_cold_region = above_region.cut_region(["obj['temperature'].in_units('K') < 3.0E5"])
below_region = data.cut_region(["(obj['cylindrical_z'].in_units('kpc') < -0.8) & np.abs(obj['cylindrical_z'].in_units('kpc') > -1.2) & (obj['cylindrical_radius'].in_units('kpc') < 2)"])
below_hot_region = below_region.cut_region(["obj['temperature'].in_units('K') > 3.0E5"])
below_cold_region = below_region.cut_region(["obj['temperature'].in_units('K') < 3.0E5"])
In [8]:
directory = '/home/aemerick/work/enzo_runs/leo_p/fiducial/sn_H2atten_H2sh/py3temp/'
data = dd.io.load(directory + 'DD1000_galaxy_data.h5')
In [9]:
print(data['time_data']['time_100'])
#print(dd.io.load(directory+'DD0200_galaxy_data.h5','/meta_data/Time'))
In [ ]:
def correct_orate(directory):
files = np.sort(glob.glob(directory + '/*galaxy_data*.h5'))
all_times = np.zeros(np.size(files))
for i,f in enumerate(files):
meta_data = dd.io.load(f,'/meta_data')
all_times[i] = meta_data['Time']
Omass = np.zeros(np.size())
times = tdata['time_100']
SNR = tdata['SNII_snr_100']
SFR = tdata['SFR_100']
Omass = np.zeros(np.size(SFR) + 1)
for i,t in enumerate(times):
index = np.argmin(np.abs(t/1.0E6-all_times))
Omass[i] = dd.io.load(files[index],
'/gas_meta_data/masses/FullBox')['O'] +\
dd.io.load(files[index],
'/gas_meta_data/masses/OutsideBox')['O']
In [29]:
def compute_loading_denominators(directory):
files = np.sort(glob.glob(directory + '/*galaxy_data*.h5'))
all_times = np.zeros(np.size(files))
#SFR = np.zeros(np.size(files))
#SNR = np.zeros(np.size(files))
#O_rate = np.zeros(np.size(files))
for i,f in enumerate(files):
meta_data = dd.io.load(f,'/meta_data')
all_times[i] = meta_data['Time']
# SFR[i] = meta_data['SFR_100']
# SNR[i] = dd.io.load(f,'/time_data/SNII_snr_100')
tdata = dd.io.load(files[-1], '/time_data')
times = tdata['time_100']
SNR = tdata['SNII_snr_100']
SFR = tdata['SFR_100']
Omass = np.zeros(np.size(SFR))
for i,t in enumerate(times):
index = np.argmin(np.abs(50.0 + t/1.0E6-all_times - all_times[0]))
Omass[i] = dd.io.load(files[index],
'/gas_meta_data/masses/FullBox')['O'] +\
dd.io.load(files[index],
'/gas_meta_data/masses/OutsideBox')['O']
Orate = np.zeros(np.size(SFR))
Orate[:-1] = (Omass[1:] - Omass[:-1]) / (times[1:] - times[:-1])
Orate[-1] = Orate[-2]
Omass = Omass / (100.0E6) # Oxygen production rate
return times, SFR, SNR, Orate, Omass
t_100, SFR_100, SNR_100, Omass_100, old_O= compute_loading_denominators(directory)
t_100 = t_100 / 1.0E6 # in Myr
f_SFR = lambda x: np.interp(x,t_100, SFR_100)
f_SNR = lambda x: np.interp(x,t_100, SNR_100)
f_Omass = lambda x: np.interp(x,t_100, Omass_100)
f_Omass_old = lambda x: np.interp(x,t_100, old_O)
print(t_100)
In [30]:
for t in datax['time']:
print(f_Omass(t), f_Omass_old(t))
In [31]:
data1 = np.genfromtxt('loading_factors.dat',names=True)
data2 = np.genfromtxt('loading_factors2.dat',names=True)
datax = {}
for k in data1.dtype.names:
datax[k] = np.array(list(data1[k]) + list(data2[k])) # ew
#print(f_Omass(datax['time']))
ofile = open("orate.dat","w")
ofile.write("#time O_rate\n")
for t in datax['time']:
ofile.write("%8.8E %8.8E\n"%(t,f_Omass(t)))
ofile.close()
In [84]:
#
# Define region
#
def compute_mass_outflow(ds, data):
dL = 200.0 * yt.units.pc
above_region = data.cut_region(["(obj['cylindrical_z'].in_units('kpc') > 0.9) & (obj['cylindrical_z'].in_units('kpc') < 1.1) & (obj['cylindrical_radius'].in_units('kpc') < 2)"])
above_hot_region = above_region.cut_region(["obj['temperature'].in_units('K') > 3.0E5"])
above_cold_region = above_region.cut_region(["obj['temperature'].in_units('K') < 3.0E5"])
above_colder_region = above_region.cut_region(["obj['temperature'].in_units('K') < 1.0E4"])
below_region = data.cut_region(["(obj['cylindrical_z'].in_units('kpc') < -0.9) & np.abs(obj['cylindrical_z'].in_units('kpc') > -1.1) & (obj['cylindrical_radius'].in_units('kpc') < 2)"])
below_hot_region = below_region.cut_region(["obj['temperature'].in_units('K') > 3.0E5"])
below_cold_region = below_region.cut_region(["obj['temperature'].in_units('K') < 3.0E5"])
below_colder_region = below_region.cut_region(["obj['temperature'].in_units('K') < 1.0E4"])
i = 0
total_vals = np.zeros(3)
metal_vals = np.zeros(3)
for above_reg, below_reg in [(above_region, below_region),
(above_hot_region,below_hot_region),
(above_cold_region,below_cold_region)]:
# (above_colder_region,below_colder_region)]:
M = above_reg['cell_mass'].to('Msun')
O = above_reg['O_Density'].value / above_reg['Density'].value * M
v = above_reg['z-velocity'].to('km/s')
select = v > 0
total_outflow_rate = np.sum(M[select]*v[select]/dL).to('Msun/yr')
metal_outflow_rate = np.sum(O[select]*v[select]/dL).to('Msun/yr')
M = below_reg['cell_mass'].to('Msun')
O = below_reg['O_Density'].value / below_reg['Density'].value * M
v = below_reg['z-velocity'].to('km/s')
select = v < 0
total_outflow_rate += np.sum(-M[select]*v[select]/dL).to('Msun/yr')
metal_outflow_rate += np.sum(-O[select]*v[select]/dL).to('Msun/yr')
total_vals[i] = total_outflow_rate
metal_vals[i] = metal_outflow_rate
i = i + 1
# now I need to compute the SFR and metal production rate
return total_vals[0], total_vals[1], total_vals[2], metal_vals[0], metal_vals[1], metal_vals[2]
def compute_energy_outflow(ds, data):
birth_mass = data['birth_mass']
birth_time = data['creation_time'].to('Myr')
pt = data['particle_type']
t = ds.current_time.to('Myr')
select = ((birth_mass > 8.0) * (birth_mass < 25.0)) * (pt > 11) * ((t - birth_time ) < 40*yt.units.Myr)
N_SN = np.size(birth_mass[select])
dL = 200 * yt.units.pc
region = None
above_region = data.cut_region(["(obj['cylindrical_z'].in_units('kpc') > 0.9) & (obj['cylindrical_z'].in_units('kpc') < 1.1) & (obj['cylindrical_radius'].in_units('kpc') < 2)"])
above_hot_region = above_region.cut_region(["obj['temperature'].in_units('K') > 3.0E5"])
above_cold_region = above_region.cut_region(["obj['temperature'].in_units('K') < 3.0E5"])
above_colder_region = above_region.cut_region(["obj['temperature'].in_units('K') < 1.0E4"])
below_region = data.cut_region(["(obj['cylindrical_z'].in_units('kpc') < -0.9) & np.abs(obj['cylindrical_z'].in_units('kpc') > -1.1) & (obj['cylindrical_radius'].in_units('kpc') < 2)"])
below_hot_region = below_region.cut_region(["obj['temperature'].in_units('K') > 3.0E5"])
below_cold_region = below_region.cut_region(["obj['temperature'].in_units('K') < 3.0E5"])
below_colder_region = below_region.cut_region(["obj['temperature'].in_units('K') < 1.0E4"])
vals = np.zeros(4)
if N_SN == 0:
print("No Supernova this time step")
return vals
i = 0
for above_reg, below_reg in [(above_region, below_region),
(above_hot_region,below_hot_region),
(above_cold_region,below_cold_region),
(above_colder_region,below_colder_region)]:
M = above_reg['cell_mass'].to('Msun')
v = above_reg['z-velocity'].to('km/s')
cV = above_reg['cell_volume'].to('cm**(3)')
KE = 0.5*M*above_reg['velocity_magnitude']**2
select = v > 0
above_E_out_therm = (v / dL * above_reg['thermal_energy'].to('erg/g') * M).to('erg/s')[select]
above_E_out_kin = (v / dL * KE).to('erg/s')[select]
M = below_reg['cell_mass'].to('Msun')
v = below_reg['z-velocity'].to('km/s')
cV = below_reg['cell_volume'].to('cm**(3)')
KE = 0.5*M*below_reg['velocity_magnitude']**2
select = v < 0
below_E_out_therm = (-v / dL * below_reg['thermal_energy'].to('erg/g') * M).to('erg/s')[select]
below_E_out_kin = (-v / dL * KE).to('erg/s')[select]
E_out_therm = np.sum(above_E_out_therm) + np.sum(below_E_out_therm)
E_out_kin = np.sum(above_E_out_kin ) + np.sum(below_E_out_kin )
E_out_tot = E_out_therm + E_out_kin
#print("Out-rate : %10.3E %10.3E %10.3E"%(E_out_therm.value, E_out_kin.value, E_out_tot.value))
#x = N_SN*yt.units.erg*1.0E51 / ((40.0*yt.units.Myr).to('s'))
x = 1.0
#print("Loading Factor: %10.3E %10.3E %10.3E"%((E_out_therm/x).value, (E_out_kin/x).value, (E_out_tot/x).value))
vals[i] = (E_out_tot/x).value * 3.154E7 # convert to erg / yr
i = i + 1
return vals
In [94]:
import glob
dsnames = np.sort(np.array(glob.glob(wdir + "DD??00/DD??00") + glob.glob(wdir+"DD??50/DD??50") +
glob.glob(wdir + "DD??25/DD??25") + glob.glob(wdir+"DD??75/DD??75")))
#print(dsnames)
dsnames = np.sort(np.array(glob.glob(wdir + "DD???5/DD???5") +glob.glob(wdir+"DD???0/DD???0")))
In [ ]:
j = 0
ds_select = dsnames[88:]
hot_loading = np.zeros(np.size(ds_select))
cold_loading = np.zeros(np.size(ds_select))
colder_loading = np.zeros(np.size(ds_select))
mass_load = np.zeros(np.size(ds_select))
mass_load_hot= np.zeros(np.size(ds_select))
mass_load_cold = np.zeros(np.size(ds_select))
metal_load = np.zeros(np.size(ds_select))
metal_load_hot = np.zeros(np.size(ds_select))
metal_load_cold = np.zeros(np.size(ds_select))
sfr = np.zeros(np.size(mass_load))
snr = np.zeros(np.size(sfr))
orate = np.zeros(np.size(snr))
loading = np.zeros(np.size(ds_select))
file = open("loading_factors_new.dat","w")
file.write("#time E_out E_hot_out E_cold_out E_colder_out M_out M_out_hot M_out_cold Metal_out Metal_out_hot Metal_out_cold SFR SNR O_rate\n")
times = np.zeros(np.size(ds_select))
for dsname in ds_select:
ds = yt.load(dsname)
data = ds.all_data()
times[j] = ds.current_time.to('Myr')
try:
x,y,z,w = compute_energy_outflow(ds,data)
except:
continue
loading[j] = x
hot_loading[j] = y
cold_loading[j] = z
colder_loading[j] = w
mass_load[j], mass_load_hot[j], mass_load_cold[j], metal_load[j], metal_load_hot[j], metal_load_cold[j] = compute_mass_outflow(ds,data)
sfr[j] = f_SFR(times[j])
snr[j] = f_SNR(times[j])
orate[j] = f_Omass(times[j])
print("%5.2f Myr: %10.3E %10.3E %10.3E %10.3E"%(times[j],loading[j],hot_loading[j],cold_loading[j],colder_loading[j]))
file.write("%5.2f %10.3E %10.3E %10.3E %10.3E %10.3E %10.3E %10.3E %10.3E %10.3E %10.3E %10.3E %10.3E %10.3E\n"%(times[j],loading[j],hot_loading[j],cold_loading[j],colder_loading[j],
mass_load[j],mass_load_hot[j], mass_load_cold[j], metal_load[j], metal_load_hot[j], metal_load_cold[j]
,sfr[j],snr[j],orate[j]))
file.flush()
j = j + 1
file.close()
In [25]:
rc('font',size=17)
fig,ax=plt.subplots()
fig.set_size_inches(8,8)
ax.plot(times,loading,lw = 3,color='black',label='Total')
ax.plot(times,hot_loading,lw=3,color='C3',label=r'(T>$3 \times 10^5$)')
ax.plot(times,cold_loading,lw=3,color='C0',label=r'(T<$3 \times 10^5$)')
ax.plot(times,colder_loading,lw=3,color='C2',label=r'(T<$1 \times 10^4$)')
select = loading>0
ax.plot(ax.get_xlim(), [np.average(hot_loading[select])]*2, lw = 2, ls = '--', color = 'C3')
ax.plot(ax.get_xlim(), [np.average(cold_loading[select])]*2, lw = 2, ls = '--', color = 'C0')
ax.plot(ax.get_xlim(), [np.average(colder_loading[select])]*2, lw = 2, ls = '--', color = 'C2')
ax.legend(loc='best',ncol=2)
ax.set_ylim(1.0E-3,15.0)
ax.set_xlim(119.0,1000.0)
ax.semilogy()
ax.set_ylabel('Energy Loading Factor')
ax.set_xlabel('Time (Myr)')
fig.savefig('energy_loading.png')
In [17]:
rc('font',size=17)
fig,ax=plt.subplots()
fig.set_size_inches(8,8)
ax.plot(times,loading,lw = 3,color='black',label='Total')
ax.plot(times,hot_loading,lw=3,color='C3',label=r'(T>$3 \times 10^5$)')
ax.plot(times,cold_loading,lw=3,color='C0',label=r'(T<$3 \times 10^5$)')
ax.plot(times,colder_loading,lw=3,color='C2',label=r'(T<$1 \times 10^4$)')
select = loading>0
ax.plot(ax.get_xlim(), [np.average(hot_loading[select])]*2, lw = 2, ls = '--', color = 'C3')
ax.plot(ax.get_xlim(), [np.average(cold_loading[select])]*2, lw = 2, ls = '--', color = 'C0')
ax.plot(ax.get_xlim(), [np.average(colder_loading[select])]*2, lw = 2, ls = '--', color = 'C2')
ax.legend(loc='best',ncol=2)
ax.set_ylim(1.0E-3,15.0)
ax.set_xlim(119.0,1000.0)
ax.semilogy()
ax.set_ylabel('Energy Loading Factor')
ax.set_xlabel('Time (Myr)')
fig.savefig('energy_loading.png')
In [88]:
data = np.genfromtxt('loading_factors.dat',names=True)
In [80]:
fig,axes = plt.subplots(1,2)
fig.set_size_inches(12,6)
ax = axes[0]
ax.plot(data['time'], data['M_out'] / data['SFR'], color = 'black', lw = 3)
ax.plot(data['time'], data['M_out_hot'] / data['SFR'], color = 'C3', lw = 3)
ax.plot(data['time'], data['M_out_cold'] / data['SFR'], color = 'C0', lw = 3)
ax.semilogy()
ax.set_xlim(119,1000.0)
ax.set_xlabel('Time (Myr)')
ax.set_ylabel('$\eta_M')
ax = axes[1]
ax.plot(data['time'], data['Metal_out'] / data['O_rate'], color = 'black', lw = 3)
ax.plot(data['time'], data['Metal_out_hot'] / data['O_rate'], color = 'C3', lw = 3)
ax.plot(data['time'], data['Metal_out_cold'] / data['O_rate'], color = 'C0', lw = 3)
ax.semilogy()
ax.set_xlim(119,1000.0)
ax.set_xlabel('Time (Myr)')
ax.set_ylabel('$\eta_Z')
ax.set_ylim(0.00001,1.0)
Out[80]:
In [83]:
fig,ax = plt.subplots()
fig.set_size_inches(6,6)
ax.plot(data['time'], data['E_out'] / data['SNR'] / 1.0E51 *3.154E7, color = 'black', lw = 3)
ax.plot(data['time'], data['E_hot_out'] / data['SNR'] *3.154E7/ 1.0E51, color = 'C3', lw = 3)
ax.plot(data['time'], data['E_cold_out'] / data['SNR']*3.154E7 / 1.0E51, color = 'C0', lw = 3)
ax.semilogy()
ax.set_xlim(119,1000.0)
ax.set_xlabel('Time (Myr)')
ax.set_ylabel('$\eta_M')
Out[83]:
In [91]:
# compute averages:
# convert to loading factors - compute averages
data = np.genfromtxt('loading_factors.dat',names=True)
loading_data = {}
for k in ['M_out','M_out_hot','M_out_cold']:
loading_data[k] = np.average( data[k] / data['SFR'])
for k in ['Metal_out','Metal_out_hot','Metal_out_cold']:
loading_data[k] = np.average( data[k] / data['O_rate'])
for k in ['E_out','E_hot_out','E_cold_out','E_colder_out']:
loading_data[k] = np.average( data[k] / (data['SNR']*1.0E51))
# compute other values
loading_data['Eta_h-Eta_c'] = loading_data['E_hot_out'] / loading_data['E_cold_out']
loading_data['e_s'] = np.average( data['E_out'] / (data['M_out']*1.989E33 )) # erg / g
loading_data['e_s_hot'] = np.average( data['E_hot_out'] / (data['M_out_cold']*1.989E33 )) # erg / g
loading_data['e_s_cold'] = np.average( data['E_cold_out'] / (data['M_out_cold']*1.989E33 )) # erg / g
loading_data['e_s_h-e_s_c'] = loading_data['e_s_hot'] / loading_data['e_s_cold']
for k in loading_data.keys():
print("%20s %8.3E"%(k,loading_data[k]))
In [93]:
for k in loading_data.keys():
print("%20s %8.3E"%(k,loading_data[k]))
In [ ]: