diff --git a/check_ics_swift.py b/check_ics_swift.py index ff3c23d..fad3198 100644 --- a/check_ics_swift.py +++ b/check_ics_swift.py @@ -8,19 +8,21 @@ Created on Tue Mar 8 11:16:06 2022 import h5py from pathlib import Path +import numpy as np -# directory = Path(r"/home/ben/sims/swiftsim/examples/agora/") +directory = Path(r"/home/ben/sims/swiftsim/examples/agora/") # directory = Path(r"/home/ben/Pictures/Gadget4/vsc4_256/") -directory = Path(r'/home/ben/sims/music-panphasia/test/') +# directory = Path(r'/home/ben/sims/music-panphasia/test/') -file = h5py.File(directory / "agora_128.hdf5", "r") -# file = h5py.File(directory / "output_0001.hdf5", "r") +# ics = h5py.File(directory / "agora_128.hdf5", "r") +file = h5py.File(directory / "output_0001.hdf5", "r") +# file_2 = h5py.File(directory / 'output_0001.hdf5', 'r') # file = h5py.File(directory / "cosmo_ics_gadget_256_pbh_vsc.hdf5", "r") -for key in file.keys(): - print(key) #for finding all header entries, which are: +# for key in file.keys(): +# print(key) #for finding all header entries, which are: Header = file["Header"] @@ -29,14 +31,64 @@ PartType1 = file["PartType1"] Units = file["Units"] # PartType2 = file['PartType2'] -print(list(Header.attrs.items())) +# print(list(Header.attrs.items())) +# print(list(ics['Header'].attrs.items())) +print(Header.attrs['Scale-factor']) -print(PartType1.keys()) -# print(PartType2.keys()) -print(PartType1['Masses'][0:10]) +# print(PartType1.keys()) +# # print(PartType2.keys()) +# print(PartType1['Masses'][0:10]) -print(PartType1['Velocities'][0:10]) +# print(PartType1['Velocities'][0:10]) + +# ic_ids = ics['PartType1']['ParticleIDs'][:] +# ic_velocities = ics['PartType1']['Velocities'][:] +# ic_absolute_vel = np.sqrt(np.sum(ic_velocities ** 2, axis=1)) +# ic_data = np.vstack([ic_ids, ic_absolute_vel]).T +# ic_positions = ics['PartType1']['Coordinates'][:] +# ic_distances = np.sqrt(np.sum(ic_positions ** 2, axis=1)) +# ic_data = np.vstack([ic_ids, ic_distances]).T + +# snap_ids = PartType1['ParticleIDs'][:] +# snap_velocities = PartType1['Velocities'][:] +# snap_absolute_vel = np.sqrt(np.sum(snap_velocities ** 2, axis=1)) +# snap_data = np.vstack([snap_ids, snap_absolute_vel]).T +# snap_data_sorted = snap_data[snap_data[:, 0].argsort()] +# snap_positions = PartType1['Coordinates'][:] +# snap_distances = np.sqrt(np.sum(snap_positions ** 2, axis=1)) +# snap_data = np.vstack([snap_ids, snap_distances]).T +# snap_data_sorted = snap_data[snap_data[:, 0].argsort()] + +# snap_2_ids = file_2['PartType1']['ParticleIDs'][:] +# snap_2_velocities = file_2['PartType1']['Velocities'][:] +# snap_2_absolute_vel = np.sqrt(np.sum(snap_2_velocities ** 2, axis=1)) +# snap_2_data = np.vstack([snap_2_ids, snap_2_absolute_vel]).T +# snap_2_data_sorted = snap_2_data[snap_2_data[:, 0].argsort()] +# snap_2_positions = file_2['PartType1']['Coordinates'][:] +# snap_2_distances = np.sqrt(np.sum(snap_2_positions ** 2, axis=1)) +# snap_2_data = np.vstack([snap_2_ids, snap_2_distances]).T +# snap_2_data_sorted = snap_2_data[snap_2_data[:, 0].argsort()] + +# velocity_difference = np.abs(ic_data[:, 1] - snap_data_sorted[:, 1]) +# mean_velocity_difference = np.mean(velocity_difference) + +# velocity_difference_2 = np.abs(snap_data_sorted[:, 1] - snap_2_data_sorted[:, 1]) +# mean_2_velocity_difference = np.mean(velocity_difference_2) + +# print(mean_velocity_difference, mean_2_velocity_difference) +# print(np.min(velocity_difference), np.max(velocity_difference)) +# print(np.min(velocity_difference_2), np.max(velocity_difference_2)) + +# distance_difference = np.abs(ic_data[:, 1] - snap_data_sorted[:, 1]) +# mean_distance_difference = np.mean(distance_difference) + +# distance_difference_2 = np.abs(snap_data_sorted[:, 1] - snap_2_data_sorted[:, 1]) +# mean_2_distance_difference = np.mean(distance_difference_2) + +# print(mean_distance_difference, mean_2_distance_difference) +# print(np.min(ic_distances), np.max(ic_distances)) +# print(np.min(snap_distances), np.max(snap_distances)) diff --git a/get_mass_zoom_halo.py b/get_mass_zoom_halo.py index 5e51b85..e7ff5a5 100644 --- a/get_mass_zoom_halo.py +++ b/get_mass_zoom_halo.py @@ -9,6 +9,10 @@ Created on Tue Mar 8 11:16:06 2022 import h5py from pathlib import Path import numpy as np +import matplotlib.pyplot as plt + +def V(r): + return 4 * np.pi * r**3 / 3 directory = Path(r"/home/ben/sims/swiftsim/examples/zoom_tests/") @@ -19,6 +23,8 @@ file = h5py.File(directory / f'output_000{snap_number}.hdf5', 'r') PartType1 = file["PartType1"] highres_groupids = PartType1['FOFGroupIDs'][:] +highres_coordinates = PartType1['Coordinates'][:] +highres_masses = PartType1['Masses'][:] total_highres_particles = len(highres_groupids) unique_groups, particle_count = np.unique(highres_groupids, return_counts=True) @@ -35,8 +41,7 @@ groups = fof_file['Groups'] # for key in groups.keys(): # print(key) - -# centres = groups['Centres'][:] +centres = groups['Centres'][:] groupids = groups['GroupIDs'][:] masses = groups['Masses'][:] number_of_members = groups['Sizes'][:] @@ -48,10 +53,61 @@ for i in range(len(groupids)-1): if np.isin(groupids[i], unique_groups): highres_members = particle_count[separate_unique_counter] contamination = (1 - highres_members / number_of_members[i]) - print(f'Group: {groupids[i]:{table_width}} | Mass: {masses[i]:{table_width + 1}.{table_width}} | Highres Members: {highres_members:{table_width}} | Total Members: {number_of_members[i]:{table_width}} | Contamination: {contamination * 100:.{table_width}}%\n') + print(f'Group: {groupids[i]:{table_width}} | Mass: {masses[i]:{table_width + 1}.{table_width}} | Highres Members: {highres_members:{table_width + 1}} | Total Members: {number_of_members[i]:{table_width}} | Contamination: {contamination * 100:.{table_width}}%\n') separate_unique_counter += 1 +main_group_origin = centres[0] +main_group_mass = masses[0] +distances = [] +# masses = [] + +for j in range(len(highres_groupids)): + if highres_groupids[j] == np.min(unique_groups): + distance = np.linalg.norm(main_group_origin - highres_coordinates[j]) + # mass = highres_masses[j] + distances.append(distance) + # masses.append(mass) + +group_radius = np.max(np.array(distances)) +normalised_distances = np.array(distances) / group_radius + +num_bins = 30 +log_radial_bins = np.geomspace(0.01, 2, num_bins) + +particle_mass = highres_masses[0] # Careful: this assumes that all highres particles have the same mass and that it's exclusively them in the main group, + # so keep an eye on its contamination and adjust accordingly! + +counts_in_radial_bins = [] +for k in range(len(log_radial_bins) - 1): + count = 0 + for member_distance in normalised_distances: + if log_radial_bins[k] <= member_distance < log_radial_bins[k + 1]: + count += 1 + volume = V(log_radial_bins[k + 1] * group_radius) - V(log_radial_bins[k] * group_radius) + count /= volume + counts_in_radial_bins.append(count) + +plot_log_radial_bins = log_radial_bins[0:len(log_radial_bins) - 1] # there is one fewer bin then there are bin borders +masses_in_radial_bins = np.array(counts_in_radial_bins) * particle_mass + +Nres = 128 +Lbox = 100 +softening = Lbox / Nres / 30 +plt.axvline(4 * softening / group_radius, linestyle='--', color='grey') + +plt.title('Density profile Auriga 6') +plt.xlabel(r'R / $\mathrm{R}_\mathrm{group}$') +plt.ylabel(r'ρ [$10^{10}\ \mathrm{M}_\odot\ /\ \mathrm{Mpc}^3$]') + +plt.loglog(plot_log_radial_bins, masses_in_radial_bins) +plt.show() + + + + + + # for i in range(10): # print(f'Group: {groupids[i]} | Mass: {masses[i]} | Size: {sizes[i]} | Centre: {centres[i]}\n') diff --git a/visualisation_swift.py b/visualisation_swift.py index 67e2d54..b1eee3f 100644 --- a/visualisation_swift.py +++ b/visualisation_swift.py @@ -41,7 +41,7 @@ def find_coordinates_to_move(minimum, maximum, ratio, x_offset, y_offset, z_offs return coordinates_to_move -directory = Path(r"/home/ben/sims/swiftsim/examples/zoom_tests/") +directory = Path(r"/home/ben/sims/swiftsim/examples/zoom_tests/auriga6_halo7_8_10/") # file = h5py.File(directory / "auriga6_halo7_8_9.hdf5", "r") @@ -54,26 +54,29 @@ for filename in sorted(directory.glob("output_*.hdf5")): highres_names = file["PartType1"]["ParticleIDs"][:] highres_velocities = file["PartType1"]["Velocities"][:] highres_masses = file['PartType1']['Masses'][:] + highres_group_ids = file['PartType1']['FOFGroupIDs'][:] highres_absolute_velo = np.sqrt(np.sum(highres_velocities ** 2, axis=1)) lowres_coordinates = file["PartType2"]["Coordinates"][:] # for cdm particles lowres_names = file["PartType2"]["ParticleIDs"][:] lowres_velocities = file["PartType2"]["Velocities"][:] lowres_masses = file['PartType2']['Masses'][:] + lowres_group_ids = file['PartType2']['FOFGroupIDs'][:] lowres_absolute_velo = np.sqrt(np.sum(lowres_velocities ** 2, axis=1)) original_coordinates = np.concatenate((highres_coordinates, lowres_coordinates)) # if "auriga" in str(filename): # original_coordinates /= 1000 - print(original_coordinates.mean()) - print(original_coordinates.min()) - print(original_coordinates.max()) + # print(original_coordinates.mean()) + # print(original_coordinates.min()) + # print(original_coordinates.max()) # print(file['Header']) # print(list(file['Units'].attrs.items())) # exit() names = np.concatenate((highres_names, lowres_names)) velocities = np.concatenate((highres_velocities, lowres_velocities)) masses = np.concatenate((highres_masses, lowres_masses)) + group_ids = np.concatenate((highres_group_ids, lowres_group_ids)) absolute_velo = np.concatenate((highres_absolute_velo, lowres_absolute_velo)) # for bla in [original_coordinates, names, velocities, masses, absolute_velo]: # print(bla.shape) @@ -86,13 +89,14 @@ for filename in sorted(directory.glob("output_*.hdf5")): velocities[::, 1], velocities[::, 2], masses, + group_ids, absolute_velo, ]).T print(original_data.shape) assert (original_coordinates == original_data[::, 0:3]).all() boundaries = Header.attrs['BoxSize'] # BoxLength for e5 boxes depends on Nres, 2.36438 for 256, 4.72876 for 512. - print(boundaries) + print(boundaries, len(highres_names)) if not boundaries.shape: boundaries = np.array([boundaries] * 3) offsets = [-1, 0, 1] @@ -191,13 +195,13 @@ for filename in sorted(directory.glob("output_*.hdf5")): # print(closest_neighbours.shape) - print(densities) - print(densities.shape) + # print(densities) + # print(densities.shape) all_data = np.column_stack([list(range(densities.shape[0])), transformed_data, densities, alt_densities]) - print(all_data.shape) - print(original_data.shape[0]) + # print(all_data.shape) + # print(original_data.shape[0]) export_data = all_data[:original_data.shape[0]] - print(export_data.shape) + # print(export_data.shape) # all_data = np.append(coordinates, velocities, axis=1) # all_data = np.column_stack((all_data, absolute_velo, names, densities)) @@ -209,5 +213,5 @@ for filename in sorted(directory.glob("output_*.hdf5")): export_data, delimiter=",", fmt="%.3f", - header="num,x,y,z,name,vx,vy,vz,masse,v,density,density_alt") + header="num,x,y,z,name,vx,vy,vz,masse,groupid,v,density,density_alt") file.close()