Skip to content

Commit 805f2ca

Browse files
glatterf42glatterf42
committed
Introduce new files for auriga evaluation and ID remapping between resolutions
1 parent 2ba0971 commit 805f2ca

5 files changed

+556
-0
lines changed

‎combine_ics.py

+298
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,298 @@
1+
#!/usr/bin/env python
2+
"""
3+
Usage:
4+
combine_ics.py input_file.0.hdf5 merged_file.hdf5 gzip_level
5+
6+
This file combines Gadget-2 type 2 (i.e. hdf5) initial condition files
7+
into a single file that can be digested by SWIFT.
8+
This has mainly be tested for DM-only (parttype1) files but also works
9+
smoothly for ICs including gas. The special case of a mass-table for
10+
the DM particles is handled. No unit conversions are applied nor are
11+
any scale-factors or h-factors changed.
12+
The script applies some compression and checksum filters to the output
13+
to save disk space.
14+
The last argument `gzip_level` is used to specify the level of compression
15+
to apply to all the fields in the file. Use 0 to cancel all coompression.
16+
The default value is `4`.
17+
18+
This file is adapted from a part of SWIFT.
19+
Copyright (C) 2016 Matthieu Schaller ([email protected])
20+
21+
All Rights Reserved.
22+
23+
This program is free software: you can redistribute it and/or modify
24+
it under the terms of the GNU Lesser General Public License as published
25+
by the Free Software Foundation, either version 3 of the License, or
26+
(at your option) any later version.
27+
28+
This program is distributed in the hope that it will be useful,
29+
but WITHOUT ANY WARRANTY; without even the implied warranty of
30+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31+
GNU General Public License for more details.
32+
33+
You should have received a copy of the GNU Lesser General Public License
34+
along with this program. If not, see <http://www.gnu.org/licenses/>.
35+
"""
36+
37+
from re import T
38+
import sys
39+
import h5py as h5
40+
import numpy as np
41+
42+
# Store the compression level
43+
gzip_level = 4
44+
if len(sys.argv) > 3:
45+
gzip_level = int(sys.argv[3])
46+
47+
# First, we need to collect some information from the master file
48+
main_file_name = str(sys.argv[1])[:-7]
49+
print("Merging snapshots files with name", main_file_name)
50+
master_file_name = main_file_name + ".0.hdf5"
51+
print("Reading master information from", master_file_name)
52+
master_file = h5.File(master_file_name, "r")
53+
grp_header = master_file["/Header"]
54+
55+
num_files = grp_header.attrs["NumFilesPerSnapshot"]
56+
tot_num_parts = grp_header.attrs["NumPart_Total"]
57+
tot_num_parts_high_word = grp_header.attrs["NumPart_Total_HighWord"]
58+
box_size = grp_header.attrs["BoxSize"]
59+
time = grp_header.attrs["Time"]
60+
hubble_param = grp_header.attrs['HubbleParam']
61+
omega0 = grp_header.attrs['Omega0']
62+
omegaLambda = grp_header.attrs['OmegaLambda']
63+
64+
# Combine the low- and high-words
65+
tot_num_parts = tot_num_parts.astype(np.int64)
66+
for i in range(6):
67+
tot_num_parts[i] += np.int64(tot_num_parts_high_word[i]) << 32
68+
69+
tot_num_parts_swift = np.copy(tot_num_parts)
70+
tot_num_parts_swift[2] += tot_num_parts_swift[3]
71+
tot_num_parts_swift[3] = 0
72+
73+
# Some basic information
74+
print("Reading", tot_num_parts, "particles from", num_files, "files.")
75+
76+
# Check whether there is a mass table
77+
DM_mass = 0.0
78+
mtable = grp_header.attrs.get("MassTable")
79+
if mtable is not None:
80+
DM_mass = grp_header.attrs["MassTable"][1] / hubble_param
81+
if DM_mass != 0.0:
82+
print("DM mass set to", DM_mass, "from the header mass table.")
83+
else:
84+
print("Reading DM mass from the particles.")
85+
86+
87+
# Create the empty file
88+
output_file_name = sys.argv[2]
89+
output_file = h5.File(output_file_name, "w-")
90+
91+
92+
# Header
93+
grp = output_file.create_group("/Header")
94+
grp.attrs["NumFilesPerSnapshot"] = 1
95+
grp.attrs["NumPart_Total"] = tot_num_parts_swift
96+
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
97+
grp.attrs["NumPart_ThisFile"] = tot_num_parts_swift
98+
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
99+
grp.attrs["BoxSize"] = box_size / hubble_param
100+
grp.attrs["Flag_Entropy_ICs"] = 0
101+
grp.attrs["Time"] = time
102+
grp.attrs['HubbleParam'] = hubble_param
103+
grp.attrs['Omega0'] = omega0
104+
grp.attrs['OmegaLambda'] = omegaLambda
105+
106+
107+
# Create the particle groups
108+
if tot_num_parts[0] > 0:
109+
grp0 = output_file.create_group("/PartType0")
110+
if tot_num_parts[1] > 0:
111+
grp1 = output_file.create_group("/PartType1")
112+
if tot_num_parts_swift[2] > 0:
113+
grp2 = output_file.create_group("/PartType2")
114+
if tot_num_parts[4] > 0:
115+
grp4 = output_file.create_group("/PartType4")
116+
if tot_num_parts[5] > 0:
117+
grp5 = output_file.create_group("/PartType5")
118+
119+
120+
# Helper function to create the datasets we need
121+
def create_set(grp, name, size, dim, dtype):
122+
if dim == 1:
123+
grp.create_dataset(
124+
name,
125+
(size,),
126+
dtype=dtype,
127+
chunks=True,
128+
compression="gzip",
129+
compression_opts=gzip_level,
130+
shuffle=True,
131+
fletcher32=True,
132+
maxshape=(size,),
133+
)
134+
else:
135+
grp.create_dataset(
136+
name,
137+
(size, dim),
138+
dtype=dtype,
139+
chunks=True,
140+
compression="gzip",
141+
compression_opts=gzip_level,
142+
shuffle=True,
143+
fletcher32=True,
144+
maxshape=(size, dim),
145+
)
146+
147+
148+
# Create the required datasets
149+
if tot_num_parts[0] > 0:
150+
create_set(grp0, "Coordinates", tot_num_parts[0], 3, "d")
151+
create_set(grp0, "Velocities", tot_num_parts[0], 3, "f")
152+
create_set(grp0, "Masses", tot_num_parts[0], 1, "f")
153+
create_set(grp0, "ParticleIDs", tot_num_parts[0], 1, "l")
154+
create_set(grp0, "InternalEnergy", tot_num_parts[0], 1, "f")
155+
create_set(grp0, "SmoothingLength", tot_num_parts[0], 1, "f")
156+
create_set(grp0, 'Potential', tot_num_parts[0], 1, 'f')
157+
158+
if tot_num_parts[1] > 0:
159+
create_set(grp1, "Coordinates", tot_num_parts[1], 3, "d")
160+
create_set(grp1, "Velocities", tot_num_parts[1], 3, "f")
161+
create_set(grp1, "Masses", tot_num_parts[1], 1, "f")
162+
create_set(grp1, "ParticleIDs", tot_num_parts[1], 1, "l")
163+
create_set(grp1, 'Potential', tot_num_parts[1], 1, 'f')
164+
165+
if tot_num_parts_swift[2] > 0:
166+
create_set(grp2, "Coordinates", tot_num_parts_swift[2], 3, "d")
167+
create_set(grp2, "Velocities", tot_num_parts_swift[2], 3, "f")
168+
create_set(grp2, "Masses", tot_num_parts_swift[2], 1, "f")
169+
create_set(grp2, "ParticleIDs", tot_num_parts_swift[2], 1, "l")
170+
create_set(grp2, 'Potential', tot_num_parts_swift[2], 1, 'f')
171+
172+
if tot_num_parts[4] > 0:
173+
create_set(grp4, "Coordinates", tot_num_parts[4], 3, "d")
174+
create_set(grp4, "Velocities", tot_num_parts[4], 3, "f")
175+
create_set(grp4, "Masses", tot_num_parts[4], 1, "f")
176+
create_set(grp4, "ParticleIDs", tot_num_parts[4], 1, "l")
177+
create_set(grp4, 'Potential', tot_num_parts[4], 1, 'f')
178+
179+
if tot_num_parts[5] > 0:
180+
create_set(grp5, "Coordinates", tot_num_parts[5], 3, "d")
181+
create_set(grp5, "Velocities", tot_num_parts[5], 3, "f")
182+
create_set(grp5, "Masses", tot_num_parts[5], 1, "f")
183+
create_set(grp5, "ParticleIDs", tot_num_parts[5], 1, "l")
184+
create_set(grp5, 'Potential', tot_num_parts[5], 1, 'f')
185+
186+
# Heavy-lifting ahead. Leave a last message.
187+
print("Datasets created in output file")
188+
189+
190+
# Special case of the non-zero mass table
191+
if DM_mass != 0.0:
192+
masses = np.ones(tot_num_parts[1], dtype=float) * DM_mass
193+
grp1["Masses"][:] = masses
194+
195+
196+
# Cumulative number of particles read/written
197+
cumul_parts = [0, 0, 0, 0, 0, 0]
198+
199+
# Loop over all the files that are part of the snapshots
200+
for f in range(num_files):
201+
202+
file_name = main_file_name + "." + str(f) + ".hdf5"
203+
file = h5.File(file_name, "r")
204+
file_header = file["/Header"]
205+
num_parts = file_header.attrs["NumPart_ThisFile"]
206+
207+
print(
208+
"Copying data from file",
209+
f,
210+
"/",
211+
num_files,
212+
": num_parts = [",
213+
num_parts[0],
214+
num_parts[1],
215+
num_parts[2],
216+
num_parts[3],
217+
num_parts[4],
218+
num_parts[5],
219+
"]",
220+
)
221+
sys.stdout.flush()
222+
223+
# Helper function to copy data
224+
def copy_grp(name_new, name_old, ptype, correct_h):
225+
full_name_new = "/PartType" + str(ptype) + "/" + name_new
226+
full_name_old = "/PartType" + str(ptype) + "/" + name_old
227+
if correct_h:
228+
output_file[full_name_new][cumul_parts[ptype] : cumul_parts[ptype] + num_parts[ptype]] = file[full_name_old] / hubble_param
229+
else:
230+
output_file[full_name_new][cumul_parts[ptype] : cumul_parts[ptype] + num_parts[ptype]] = file[full_name_old]
231+
232+
def copy_grp_same_name(name, ptype, correct_h):
233+
copy_grp(name, name, ptype, correct_h)
234+
235+
def copy_grp_pt3(name, correct_h):
236+
full_name_new = "/PartType2/" + name
237+
full_name_old = "/PartType3/" + name
238+
if correct_h:
239+
output_file[full_name_new][cumul_parts[2] : cumul_parts[2] + num_parts[3]] = file[full_name_old] / hubble_param
240+
else:
241+
output_file[full_name_new][cumul_parts[2] : cumul_parts[2] + num_parts[3]] = file[full_name_old]
242+
243+
if num_parts[0] > 0:
244+
copy_grp_same_name("Coordinates", 0, True)
245+
copy_grp_same_name("Velocities", 0, False)
246+
copy_grp_same_name("Masses", 0, False)
247+
copy_grp_same_name("ParticleIDs", 0, False)
248+
copy_grp_same_name("InternalEnergy", 0, False)
249+
copy_grp_same_name("SmoothingLength", 0, False)
250+
copy_grp_same_name('Potential', 0, False) ######## Careful: I don't actually know the units of the Potential, so I don't know if I should correct it.
251+
252+
if num_parts[1] > 0:
253+
copy_grp_same_name("Coordinates", 1, True)
254+
copy_grp_same_name("Velocities", 1, False)
255+
copy_grp_same_name("ParticleIDs", 1, False)
256+
if DM_mass == 0.0: # Do not overwrite values if there was a mass table
257+
copy_grp_same_name("Masses", 1, False)
258+
copy_grp_same_name('Potential', 1, False)
259+
260+
if num_parts[2] > 0:
261+
copy_grp_same_name("Coordinates", 2, True)
262+
copy_grp_same_name("Velocities", 2, False)
263+
copy_grp_same_name("ParticleIDs", 2, False)
264+
copy_grp_same_name("Masses", 2, False)
265+
copy_grp_same_name('Potential', 2, False)
266+
267+
# Need to update part counter for pt2 already here, so we append 3 in correct place
268+
cumul_parts[2] += num_parts[2]
269+
270+
if num_parts[3] > 0:
271+
copy_grp_pt3("Coordinates", True)
272+
copy_grp_pt3("Velocities", False)
273+
copy_grp_pt3("ParticleIDs", False)
274+
copy_grp_pt3("Masses", False)
275+
copy_grp_pt3('Potential', False)
276+
277+
if num_parts[4] > 0:
278+
copy_grp_same_name("Coordinates", 4, True)
279+
copy_grp_same_name("Velocities", 4, False)
280+
copy_grp_same_name("Masses", 4, False)
281+
copy_grp_same_name("ParticleIDs", 4, False)
282+
copy_grp_same_name('Potential', 4, False)
283+
284+
if num_parts[5] > 0:
285+
copy_grp_same_name("Coordinates", 5, True)
286+
copy_grp_same_name("Velocities", 5, False)
287+
copy_grp_same_name("Masses", 5, False)
288+
copy_grp_same_name("ParticleIDs", 5, False)
289+
copy_grp_same_name('Potential', 5, False)
290+
291+
cumul_parts[0] += num_parts[0]
292+
cumul_parts[1] += num_parts[1]
293+
cumul_parts[2] += num_parts[3] # Need to adjust for added pt-3
294+
cumul_parts[4] += num_parts[4]
295+
cumul_parts[5] += num_parts[5]
296+
file.close()
297+
298+
print("All done! SWIFT is waiting.")

‎hist2d_auriga6.py

+69
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
#!/usr/bin/env python3
2+
# -*- coding: utf-8 -*-
3+
"""
4+
Created on Wed Apr 27 11:55:06 2022
5+
6+
@author: ben
7+
"""
8+
9+
import h5py
10+
from pathlib import Path
11+
import numpy as np
12+
import matplotlib.pyplot as plt
13+
from matplotlib.colors import LogNorm
14+
15+
directory = Path(r"/home/ben/sims/swiftsim/examples/zoom_tests/auriga6_halo7_8_10")
16+
Nres = 128 # 233.45 #for arj
17+
Lbox = 100
18+
19+
snap_number = 7 # 7 for our tests, 0 for arj
20+
21+
fof_file = h5py.File(directory / f"fof_output_000{snap_number}.hdf5", "r")
22+
file = h5py.File(directory / f'output_000{snap_number}.hdf5', 'r')
23+
24+
PartType1 = file["PartType1"]
25+
highres_groupids = PartType1['FOFGroupIDs'][:]
26+
highres_coordinates = np.array(PartType1['Coordinates'][:])
27+
highres_masses = PartType1['Masses'][:]
28+
29+
total_highres_particles = len(highres_groupids)
30+
unique_groups, particle_count = np.unique(highres_groupids, return_counts=True)
31+
# print(f'Group IDs: {unique_groups}, Number in them: {particle_count}')
32+
33+
# all particles belonging to no group have group id 2147483647
34+
35+
highres_x, highres_y, highres_z = highres_coordinates[:, 0], highres_coordinates[:, 1], highres_coordinates[:, 2]
36+
37+
groups = fof_file['Groups']
38+
39+
# for key in groups.keys():
40+
# print(key)
41+
42+
centres = groups['Centres'][:]
43+
groupids = groups['GroupIDs'][:]
44+
masses = groups['Masses'][:]
45+
number_of_members = groups['Sizes'][:]
46+
47+
# table_width = 4
48+
49+
# separate_unique_counter = 0
50+
# # for i in range(len(groupids)-1):
51+
# for i in range(11):
52+
# if np.isin(groupids[i], unique_groups):
53+
# highres_members = particle_count[separate_unique_counter]
54+
# contamination = (1 - highres_members / number_of_members[i])
55+
# 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')
56+
# separate_unique_counter += 1
57+
58+
main_group_origin = centres[np.min(unique_groups) - 1]
59+
main_group_mass = masses[np.min(unique_groups) - 1]
60+
origin_x, origin_y, origin_z = main_group_origin[0], main_group_origin[1], main_group_origin[2]
61+
print(np.mean(highres_masses))
62+
63+
plt.title('Histogram of Auriga 6 Centre Levelmax 10')
64+
plt.xlabel('x [Mpc]')
65+
plt.ylabel('y [Mpc]')
66+
i = np.logical_and( highres_z>origin_z-0.5, highres_z<origin_z+0.5 )
67+
plt.hist2d(x=highres_x[i], y=highres_y[i], bins=256, range=[[origin_x - 0.5, origin_x + 0.5], [origin_y - 0.5, origin_y + 0.5]], weights=highres_masses[i], norm=LogNorm())
68+
plt.colorbar()
69+
plt.show()

0 commit comments

Comments
 (0)