From 3c4db1658e20825f50e0815d1ce8717e83d7ec00 Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Fri, 4 Feb 2022 16:38:14 +0100 Subject: [PATCH] Add clump stuff + solve issue with pdf --- snapshotprocessor.py | 163 +++++++++++++++++++++++++++++-------------- 1 file changed, 112 insertions(+), 51 deletions(-) diff --git a/snapshotprocessor.py b/snapshotprocessor.py index 2ce2a36..5ffb6fc 100644 --- a/snapshotprocessor.py +++ b/snapshotprocessor.py @@ -5,6 +5,7 @@ import tables import pickle import astropy.units as u import pandas as pd +from struct import pack from skimage.morphology import medial_axis @@ -28,6 +29,7 @@ from pymses.analysis import ( splatting, ) from pymses.filters import CellsToPoints, RegionFilter +from pymses.sources.hop.hop import HOP from fil_finder import FilFinder2D from scipy import fft @@ -634,7 +636,7 @@ class SnapshotProcessor(HDF5Container): """ Azimuthal velocity """ return self.oct_getter_vect_phi(dset, "vel") - def _slice(self, getter, ax_los="z", z=0.0, unit=U.none): + def slice(self, getter, ax_los="z", z=0.0, unit=U.none): """ Slice process function. Return a slice of the source box. @@ -662,7 +664,7 @@ class SnapshotProcessor(HDF5Container): datamap = slicing.SliceMap(self._amr, self._cam[ax_los], op, z=z) return datamap.map.T - def _ax_avg(self, getter, ax_los, unit=U.none, mass_weighted=True, surf_qty=False): + def ax_avg(self, getter, ax_los, unit=U.none, mass_weighted=True, surf_qty=False): """ Map of the average of a quantity (given by getter) along an axis (ax_los) Returns 2D array if getter returns a scalar quantity @@ -741,7 +743,7 @@ class SnapshotProcessor(HDF5Container): return df.groupby("axis").mean().values[:, 0] - def _sum(self, getter, mass_weighted=True): + def sum(self, getter, mass_weighted=True): """ Global sum of the quantity returned by getter (variable must be extensive) Returns a scalar (or a vector if the quantity returned by getter is a getter, eg. speed) @@ -759,7 +761,7 @@ class SnapshotProcessor(HDF5Container): self.unload_cells() return data - def _vol_avg(self, getter, mass_weighted=True): + def vol_avg(self, getter, mass_weighted=True): """ Global volumic (or mass_weighted) average of the quantity returned by getter Returns a scalar (or a vector if the quantity returned by getter is a getter, eg. speed) @@ -777,7 +779,7 @@ class SnapshotProcessor(HDF5Container): self.unload_cells() return data - def _vol_pdf(self, getter, bins=100, logbins=False, weight_func=vol_func): + def vol_pdf(self, getter, bins=100, logbins=False, weight_func=vol_func): self.load_cells() data = getter(self.cells) if logbins: @@ -786,7 +788,7 @@ class SnapshotProcessor(HDF5Container): if self.params.process.unload_cells: self.unload_cells() - values, edges = np.histogram(data, bins, weights=weights) + values, edges = np.histogram(data, bins, weights=weights, density=True) centers = 0.5 * (edges[1:] + edges[:-1]) return (np.stack([values, centers]), {"logbins": logbins}) @@ -804,7 +806,7 @@ class SnapshotProcessor(HDF5Container): self.unload_cells() return ({"rho": centers, "B": B_mean}, {"logbins": logbins}) - def _mean_by_bins( + def mean_by_bins( self, name_x, name_y, ax_los, group="/maps/", bins=100, logbins=True ): """ @@ -885,7 +887,7 @@ class SnapshotProcessor(HDF5Container): dot_prod = np.einsum("ij,ij->i", vel_fluct, dset["Br"]) return np.abs(dot_prod) / (v_norm * B_norm) - return self._vol_pdf(getter_cos_vfluct_B) + return self.vol_pdf(getter_cos_vfluct_B) def _mwa_sigma(self, axes=["x", "y", "z"]): mw_speed = self.get_value("/globals/mwa_speed") @@ -905,7 +907,7 @@ class SnapshotProcessor(HDF5Container): sigma_squared = sigma_squared + sigma_sq_ax return sigma_squared - return np.sqrt(self._vol_avg(getter, mass_weighted=True)) + return np.sqrt(self.vol_avg(getter, mass_weighted=True)) def _coldens(self, ax_los): datamap = self._rt.process(self._cam[ax_los], surf_qty=True) @@ -1171,9 +1173,9 @@ class SnapshotProcessor(HDF5Container): dmap = self.get_value("/maps/" + name + "_" + ax_los) if mass_weighted: - avg_map = self.get_value("/maps/mwavg_map_" + name + "_" + ax_los) + avg_map = self.get_value("/maps/rad_mwavg_map_" + name + "_" + ax_los) else: - avg_map = self.get_value("/maps/avg_map_" + name + "_" + ax_los) + avg_map = self.get_value("/maps/rad_avg_map_" + name + "_" + ax_los) return dmap / avg_map def _dispersion_rad(self, name, ax_los="z"): @@ -1189,7 +1191,7 @@ class SnapshotProcessor(HDF5Container): # 1. Get full storage names for the map and its azimuthal avererage map_name = f"/maps/{name}_{ax_los}" - avg_map_name = f"/maps/avg_map_{name}_{ax_los}" + avg_map_name = f"/maps/rad_avg_map_{name}_{ax_los}" # 2. Get maps from the storage dmap = self.get_value(map_name) @@ -1236,7 +1238,7 @@ class SnapshotProcessor(HDF5Container): return np.stack([values, centers]) def _fit_pdf(self, name, ax_los="z"): - pdf = self.get_value("/hist/pdf_" + name + "_" + ax_los) + pdf = self.get_value("/hist/fluct_pdf_" + name + "_" + ax_los) values, centers = pdf mask_fit = ( (centers > self.params.pdf.xmin_fit) @@ -1260,7 +1262,7 @@ class SnapshotProcessor(HDF5Container): # Mean part - T_avg = self.get_value("/maps/avg_map_T_mwavg_z") + T_avg = self.get_value("/maps/rad_avg_map_T_mwavg_z") radial_bins = self.get_value("/radial/radial_bins_" + ax_los) mean_bin_vr = self.get_value("/radial/rad_avg_" + "slice_velr" + "_" + ax_los) @@ -1299,7 +1301,7 @@ class SnapshotProcessor(HDF5Container): return alpha alpha_f = ( - self._ax_avg(getter_alpha_num, "z", unit=U.none, mass_weighted=True) / T_avg + self.ax_avg(getter_alpha_num, "z", unit=U.none, mass_weighted=True) / T_avg ) # alpha @@ -1310,8 +1312,8 @@ class SnapshotProcessor(HDF5Container): "Map of the gravitational contribution to the Shakura&Sunaev alpha parameter for disks" assert ax_los == "z" - T_avg = self.get_value("/maps/avg_map_T_mwavg_z") - coldens = self.get_value("/maps/avg_map_coldens_z") + T_avg = self.get_value("/maps/rad_avg_map_T_mwavg_z") + coldens = self.get_value("/maps/rad_avg_map_coldens_z") def getter_alpha_grav(dset): r2 = np.sum((self.lbox * self.oct_getter_pos_disk(dset)) ** 2, axis=2) @@ -1321,7 +1323,7 @@ class SnapshotProcessor(HDF5Container): gphi = self.oct_getter_vect_phi(dset, "g") return gr * gphi / (4 * np.pi * self.G) - alpha_g = self._ax_avg(getter_alpha_grav, "z", unit=U.none, surf_qty=True) / ( + alpha_g = self.ax_avg(getter_alpha_grav, "z", unit=U.none, surf_qty=True) / ( coldens * T_avg ) @@ -1515,6 +1517,70 @@ class SnapshotProcessor(HDF5Container): return {"gfil": gfil, "Rfil": Rfil, "fPfil": fPfil, "dvr": dvr_fil} + def make_clump_hop(self, threshold_density=1e2): + """ + Apply HOP algorithm + + Args: + threshold_density (float): select only cells over threshold + """ + + # Selection of cells + mask = self.cells["rho"] * self.info["unit_density"].express(U.H_cc) > threshold_density + ncells = np.sum(mask) + + # fill the matrice with ID, x,y,z and masses of particles + cells_group = np.zeros((ncells, 6)) + cells_group[:,0] = np.arange(ncells) # index + position = self.cells["pos"][mask] * self.info["unit_length"].express(U.pc) + cells_group[:,1:4] = position # position + density = self.cells["rho"][mask] * self.info["unit_density"].express(U.Msun / U.pc**3) + size = self.cells["dx"][mask]*self.info["unit_length"].express(U.pc) + cells_group[:,4] = density * size**3 # mass + + # save file.txt + head = str(ncells) + np.savetxt( + self.filename[:-3] + '_hop.txt', + cells_group[:, :-1], + fmt='%10d %.10e %.10e %.10e %.10e', + header=head, + delimiter=' ', + comments=' ' + ) + + # save file.den + f = open(self.filename[:-3] + '_hop.den','wb') + f.write(pack('I', ncells)) + self.cells["rho"][mask].astype("f").tofile(f) + f.close() + + # exec HOP algo + h = HOP(self.filename[:-3] + '_hop.txt', os.path.dirname(self.filename)) + h.process_hop() + + # get the igroup array + group_ids = h.get_group_ids() + + # sort it and apply the sorting to the coordinates + # this means that the particules of group 1 are written first then of group 2 etc... + ind_sort = np.argsort(group_ids) + cells_group = val_mat[ind_sort] + cells_group[6] = group_ids[ind_sort] + + # Make it a pandas' DataFrame + cells_group = pd.DataFrame(cells_group, header=["id", "x", "y", "z", "mass", "group"]) + + self.clumps = cells_group + + return cells_group + + def clump_properties(self, nograv=False, log_min=None, log_max=None): + + cells_group = self.make_clump_hop() + cells_group.groupby("group") + + def def_rules(self): self.rules = { @@ -1529,7 +1595,7 @@ class SnapshotProcessor(HDF5Container): "T_mwavg": Rule( self, partial( - self._ax_avg, + self.ax_avg, getter_T, mass_weighted=True, unit=self.info["unit_pressure"] / self.info["unit_density"], @@ -1702,7 +1768,7 @@ class SnapshotProcessor(HDF5Container): # PDF "rho_pdf": Rule( self, - partial(self._vol_pdf, partial(simple_getter, "rho"), logbins=True), + partial(self.vol_pdf, partial(simple_getter, "rho"), logbins=True), "Global rho-PDF", "/hist", unit=self.info["unit_density"], @@ -1710,7 +1776,7 @@ class SnapshotProcessor(HDF5Container): "rho_pdf_mw": Rule( self, partial( - self._vol_pdf, + self.vol_pdf, partial(simple_getter, "rho"), weight_func=mass_func, logbins=True, @@ -1721,18 +1787,11 @@ class SnapshotProcessor(HDF5Container): ), "T_pdf": Rule( self, - partial(self._vol_pdf, getter_T, logbins=True), + partial(self.vol_pdf, getter_T, logbins=True), "Global T-PDF", "/hist", unit=self.info["unit_pressure"] / self.info["unit_density"], ), - "P_pdf": Rule( - self, - partial(self._vol_pdf, getter_T, logbins=True), - "Global P-PDF", - "/hist", - unit=self.info["unit_pressure"], - ), "cos_pdf": Rule( self, partial(self.cos_vfluct_B), @@ -1782,14 +1841,14 @@ class SnapshotProcessor(HDF5Container): ), "mass": Rule( self, - partial(self._sum, mass_func, mass_weighted=False), + partial(self.sum, mass_func, mass_weighted=False), "Total mass", "/globals", unit=self.info["unit_density"] * self.info["unit_length"] ** 3, ), "mwa_speed": Rule( self, - partial(self._vol_avg, partial(simple_getter, "vel")), + partial(self.vol_avg, partial(simple_getter, "vel")), "Mass weighted speed average", "/globals", unit=self.info["unit_velocity"], @@ -1804,14 +1863,14 @@ class SnapshotProcessor(HDF5Container): ), "mwa_B_int": Rule( self, - partial(self._vol_avg, getter_B_int), + partial(self.vol_avg, getter_B_int), "Mass weighted Magnetic intensity average", "/globals", unit=self.info["unit_mag"], ), } - averageables = [ + averageable_maps = [ "coldens", "Q", "T", @@ -1829,7 +1888,7 @@ class SnapshotProcessor(HDF5Container): self.rules["slice_" + name] = Rule( self, - partial(self._slice, getter, z=0.0, unit=unit), + partial(self.slice, getter, z=0.0, unit=unit), "{} slice".format(name), "/maps", unit=unit, @@ -1837,7 +1896,7 @@ class SnapshotProcessor(HDF5Container): self.rules[name + "_mwavg"] = Rule( self, - partial(self._ax_avg, oct_getter, mass_weighted=True, unit=unit), + partial(self.ax_avg, oct_getter, mass_weighted=True, unit=unit), "Ax mass-weighted averaged {}".format(name), "/maps", unit=unit, @@ -1845,15 +1904,15 @@ class SnapshotProcessor(HDF5Container): self.rules[name + "_avg"] = Rule( self, - partial(self._ax_avg, oct_getter, mass_weighted=False, unit=unit), + partial(self.ax_avg, oct_getter, mass_weighted=False, unit=unit), "Ax averaged {}".format(name), "/maps", unit=unit, ) - averageables.append("slice_" + name) - averageables.append(name + "_mwavg") - averageables.append(name + "_avg") + averageable_maps.append("slice_" + name) + averageable_maps.append(name + "_mwavg") + averageable_maps.append(name + "_avg") # special for vectors if field in ["g", "vel", "Bl", "Br"]: @@ -1891,10 +1950,12 @@ class SnapshotProcessor(HDF5Container): generic_rule(field, partial(simple_getter, field), self.unit_key[field]) # radial average and other - for name in averageables: + for name in averageable_maps: unit = self.rules[name].unit + + self.rules["rad_avg_" + name] = Rule( self, partial(self._rad_avg, name), @@ -1911,7 +1972,7 @@ class SnapshotProcessor(HDF5Container): dependencies=["coldens", "radial_centers", "bins_on_map", name], unit=unit, ) - self.rules["avg_map_" + name] = Rule( + self.rules["rad_avg_map_" + name] = Rule( self, partial(self._rad_avg_map, name), "Interpolated map of azimuthal average of {}".format(name), @@ -1919,7 +1980,7 @@ class SnapshotProcessor(HDF5Container): dependencies=["radial_centers", "bins_on_map", "rad_avg_" + name], unit=unit, ) - self.rules["mwavg_map_" + name] = Rule( + self.rules["rad_mwavg_map_" + name] = Rule( self, partial(self._rad_avg_map, name, mass_weighted=True), "Interpolated map of azimuthal average of {}".format(name), @@ -1932,14 +1993,14 @@ class SnapshotProcessor(HDF5Container): partial(self._fluct_map, name), "Fluctuation wrt to average of {}".format(name), "/maps", - dependencies=[name, "avg_map_" + name], + dependencies=[name, "rad_avg_map_" + name], ) self.rules["dispersion_rad_" + name] = Rule( self, partial(self._dispersion_rad, name), "radial RMS of {}".format(name), "/radial", - dependencies=[name, "avg_map_" + name], + dependencies=[name, "rad_avg_map_" + name], unit=unit, ) self.rules["mwfluct_" + name] = Rule( @@ -1947,16 +2008,16 @@ class SnapshotProcessor(HDF5Container): partial(self._fluct_map, name, mass_weigthed=True), "Fluctuation wrt to mass weighted average of {}".format(name), "/maps", - dependencies=[name, "mwavg_map_" + name], + dependencies=[name, "rad_mwavg_map_" + name], ) - self.rules["pdf_" + name] = Rule( + self.rules["fluct_pdf_" + name] = Rule( self, partial(self._rad_fluct_pdf, name), "Probability density function of {} fluctuations".format(name), "/hist", dependencies=["rr", "fluct_" + name], ) - self.rules["mwpdf_" + name] = Rule( + self.rules["fluct_mwpdf_" + name] = Rule( self, partial(self._rad_fluct_pdf, name, mass_weigthed=True), "Probability density function of {} mass weighted fluctuations".format( @@ -1965,19 +2026,19 @@ class SnapshotProcessor(HDF5Container): "/hist", dependencies=["rr", "mwfluct_" + name], ) - self.rules["fit_pdf_" + name] = Rule( + self.rules["fit_fluct_pdf_" + name] = Rule( self, partial(self._fit_pdf, name), "Fit the PDF of {} fluctuations".format(name), "/hist", - dependencies=["pdf_" + name], + dependencies=["fluct_pdf_" + name], ) - for name_bin in averageables: + for name_bin in averageable_maps: unit_bin = self.rules[name_bin].unit if name_bin is not name: self.rules["mbb_" + name + "_" + name_bin] = Rule( self, - partial(self._mean_by_bins, name_bin, name), + partial(self.mean_by_bins, name_bin, name), "Mean of {} by bins of {}".format(name, name_bin), "/hist", dependencies=[name, name_bin],