Commit 7cd656bd by Stefan Arridge

### Added READMEs to 'cooling halo' and 'cooling halo with spin'.

parent f83f665c
 ; ;To make the initial conditions we distribute gas particles randomly in a cube ;with a side length twice that of the virial radius. The density profile ;of the gas is proportional to r^(-2) where r is the distance from the centre ;of the cube. ; ;The parameter v_rot (in makeIC.py and cooling.yml) sets the circular velocity ;of the halo, and by extension, the viral radius, viral mass, and the ;internal energy of the gas such that hydrostatic equilibrium is achieved. ; ;While the system is initially in hydrostatic equilibrium, the cooling of the ;gas means that the halo will collapse. ; ;To run this example, make such that the code is compiled with either the ;isothermal potential or softened isothermal potential, and 'const_lambda' ;cooling, set in src/const.h. In ;the latter case, a (small) value of epsilon needs to be set in cooling.yml. ;0.1 kpc should work well. ; ;The plotting scripts produce a plot of the density, internal energy and radial ;velocity profile for each snapshot. test_energy_conservation.py shows the ;evolution of energy with time. These can be used to check if the example ;has run properly. ; ; \ No newline at end of file
 import numpy as np import h5py as h5 import matplotlib.pyplot as plt import sys def do_binning(x,y,x_bin_edges): #x and y are arrays, where y = f(x) #returns number of elements of x in each bin, and the total of the y elements corresponding to those x values n_bins = x_bin_edges.size - 1 count = np.zeros(n_bins) y_totals = np.zeros(n_bins) for i in range(n_bins): ind = np.intersect1d(np.where(x > bin_edges[i])[0],np.where(x <= bin_edges[i+1])[0]) count[i] = ind.size binned_y = y[ind] y_totals[i] = np.sum(binned_y) return(count,y_totals) #for the plotting max_r = float(sys.argv[1]) n_radial_bins = int(sys.argv[2]) n_snaps = int(sys.argv[3]) #some constants OMEGA = 0.3 # Cosmological matter fraction at z = 0 PARSEC_IN_CGS = 3.0856776e18 KM_PER_SEC_IN_CGS = 1.0e5 CONST_G_CGS = 6.672e-8 CONST_m_H_CGS = 1.67e-24 h = 0.67777 # hubble parameter gamma = 5./3. eta = 1.2349 H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS) #read some header/parameter information from the first snapshot filename = "CoolingHalo_000.hdf5" f = h5.File(filename,'r') params = f["Parameters"] unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"]) unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"]) unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"]) unit_time_cgs = unit_length_cgs / unit_velocity_cgs v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"]) v_c_cgs = v_c * unit_velocity_cgs header = f["Header"] N = header.attrs["NumPart_Total"][0] box_centre = np.array(header.attrs["BoxSize"]) #calculate r_vir and M_vir from v_c r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA)) M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS for i in range(n_snaps): filename = "CoolingHalo_%03d.hdf5" %i f = h5.File(filename,'r') coords_dset = f["PartType0/Coordinates"] coords = np.array(coords_dset) #translate coords by centre of box header = f["Header"] snap_time = header.attrs["Time"] snap_time_cgs = snap_time * unit_time_cgs coords[:,0] -= box_centre[0]/2. coords[:,1] -= box_centre[1]/2. coords[:,2] -= box_centre[2]/2. radius = np.sqrt(coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2) radius_cgs = radius*unit_length_cgs radius_over_virial_radius = radius_cgs / r_vir_cgs #get the internal energies vel_dset = f["PartType0/Velocities"] vel = np.array(vel_dset) #make dimensionless vel /= v_c r = radius_over_virial_radius #find radial component of velocity v_r = np.zeros(r.size) for j in range(r.size): v_r[j] = -np.dot(coords[j,:],vel[j,:])/radius[j] bin_edges = np.linspace(0,max_r,n_radial_bins + 1) (hist,v_r_totals) = do_binning(r,v_r,bin_edges) bin_widths = bin_edges[1] - bin_edges[0] radial_bin_mids = np.linspace(bin_widths / 2. , max_r - bin_widths / 2. , n_radial_bins) binned_v_r = v_r_totals / hist #calculate cooling radius #r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs) plt.plot(radial_bin_mids,binned_v_r,'ko',label = "Average radial velocity in shell") #plt.plot((0,1),(1,1),label = "Analytic Solution") #plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,2),'r',label = "Cooling radius") plt.legend(loc = "upper right") plt.xlabel(r"$r / r_{vir}$") plt.ylabel(r"$v_r / v_c$") plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c)) plt.ylim((0,2)) plot_filename = "./plots/radial_velocity_profile/velocity_profile_%03d.png" %i plt.savefig(plot_filename,format = "png") plt.close()
 ; ;To make the initial conditions we distribute gas particles randomly in a cube ;with a side length twice that of the virial radius. The density profile ;of the gas is proportional to r^(-2) where r is the distance from the centre ;of the cube. ; ;The parameter v_rot (in makeIC.py and cooling.yml) sets the circular velocity ;of the halo, and by extension, the viral radius, viral mass, and the ;internal energy of the gas such that hydrostatic equilibrium is achieved. ; ;The gas is given some angular momentum about the z-axis. This is defined by ;the 'spin_lambda' variable in makeIC.py ; ;While the system is initially in hydrostatic equilibrium, the cooling of the ;gas and the non-zero angular momentum means that the halo will collapse into ;a spinning disc. ; ;To run this example, make such that the code is compiled with either the ;isothermal potential or softened isothermal potential, and 'const_lambda' ;cooling, set in src/const.h. In ;the latter case, a (small) value of epsilon needs to be set in cooling.yml. ;0.1 kpc should work well. ; ;The plotting scripts produce a plot of the density, internal energy and radial ;velocity profile for each snapshot. test_energy_conservation.py shows the ;evolution of energy with time. These can be used to check if the example ;has run properly. ; ; \ No newline at end of file
 import numpy as np import h5py as h5 import matplotlib.pyplot as plt import sys def do_binning(x,y,x_bin_edges): #x and y are arrays, where y = f(x) #returns number of elements of x in each bin, and the total of the y elements corresponding to those x values n_bins = x_bin_edges.size - 1 count = np.zeros(n_bins) y_totals = np.zeros(n_bins) for i in range(n_bins): ind = np.intersect1d(np.where(x > bin_edges[i])[0],np.where(x <= bin_edges[i+1])[0]) count[i] = ind.size binned_y = y[ind] y_totals[i] = np.sum(binned_y) return(count,y_totals) #for the plotting max_r = float(sys.argv[1]) n_radial_bins = int(sys.argv[2]) n_snaps = int(sys.argv[3]) #some constants OMEGA = 0.3 # Cosmological matter fraction at z = 0 PARSEC_IN_CGS = 3.0856776e18 KM_PER_SEC_IN_CGS = 1.0e5 CONST_G_CGS = 6.672e-8 CONST_m_H_CGS = 1.67e-24 h = 0.67777 # hubble parameter gamma = 5./3. eta = 1.2349 H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS) #read some header/parameter information from the first snapshot filename = "CoolingHalo_000.hdf5" f = h5.File(filename,'r') params = f["Parameters"] unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"]) unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"]) unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"]) unit_time_cgs = unit_length_cgs / unit_velocity_cgs v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"]) v_c_cgs = v_c * unit_velocity_cgs header = f["Header"] N = header.attrs["NumPart_Total"][0] box_centre = np.array(header.attrs["BoxSize"]) #calculate r_vir and M_vir from v_c r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA)) M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS for i in range(n_snaps): filename = "CoolingHalo_%03d.hdf5" %i f = h5.File(filename,'r') coords_dset = f["PartType0/Coordinates"] coords = np.array(coords_dset) #translate coords by centre of box header = f["Header"] snap_time = header.attrs["Time"] snap_time_cgs = snap_time * unit_time_cgs coords[:,0] -= box_centre[0]/2. coords[:,1] -= box_centre[1]/2. coords[:,2] -= box_centre[2]/2. radius = np.sqrt(coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2) radius_cgs = radius*unit_length_cgs radius_over_virial_radius = radius_cgs / r_vir_cgs #get the internal energies vel_dset = f["PartType0/Velocities"] vel = np.array(vel_dset) #make dimensionless vel /= v_c r = radius_over_virial_radius #find radial component of velocity v_r = np.zeros(r.size) for j in range(r.size): v_r[j] = -np.dot(coords[j,:],vel[j,:])/radius[j] bin_edges = np.linspace(0,max_r,n_radial_bins + 1) (hist,v_r_totals) = do_binning(r,v_r,bin_edges) bin_widths = bin_edges[1] - bin_edges[0] radial_bin_mids = np.linspace(bin_widths / 2. , max_r - bin_widths / 2. , n_radial_bins) binned_v_r = v_r_totals / hist #calculate cooling radius #r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs) plt.plot(radial_bin_mids,binned_v_r,'ko',label = "Average radial velocity in shell") #plt.plot((0,1),(1,1),label = "Analytic Solution") #plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,2),'r',label = "Cooling radius") plt.legend(loc = "upper right") plt.xlabel(r"$r / r_{vir}$") plt.ylabel(r"$v_r / v_c$") plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c)) plt.ylim((0,2)) plot_filename = "./plots/radial_velocity_profile/velocity_profile_%03d.png" %i plt.savefig(plot_filename,format = "png") plt.close()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment