noe/pipeline
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Add clump stuff + solve issue with pdf

Browse Source
This commit is contained in:
Noe Brucy
2022-02-04 16:38:14 +01:00
parent fd367b7a28
commit 3c4db1658e
1 changed files with 112 additions and 51 deletions
Download Patch File Download Diff File Expand all files Collapse all files
snapshotprocessor.py
+112 -51
View File
@@ -5,6 +5,7 @@ import tables
import pickle import pickle
import astropy.units as u import astropy.units as u
import pandas as pd import pandas as pd
from struct import pack
from skimage.morphology import medial_axis from skimage.morphology import medial_axis
@@ -28,6 +29,7 @@ from pymses.analysis import (
splatting, splatting,
) )
from pymses.filters import CellsToPoints, RegionFilter from pymses.filters import CellsToPoints, RegionFilter
from pymses.sources.hop.hop import HOP
from fil_finder import FilFinder2D from fil_finder import FilFinder2D
from scipy import fft from scipy import fft
@@ -634,7 +636,7 @@ class SnapshotProcessor(HDF5Container):
""" Azimuthal velocity """ """ Azimuthal velocity """
return self.oct_getter_vect_phi(dset, "vel") 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. Slice process function.
Return a slice of the source box. 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) datamap = slicing.SliceMap(self._amr, self._cam[ax_los], op, z=z)
return datamap.map.T 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) Map of the average of a quantity (given by getter) along an axis (ax_los)
Returns 2D array if getter returns a scalar quantity Returns 2D array if getter returns a scalar quantity
@@ -741,7 +743,7 @@ class SnapshotProcessor(HDF5Container):
return df.groupby("axis").mean().values[:, 0] 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) 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) 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() self.unload_cells()
return data 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 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) 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() self.unload_cells()
return data 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() self.load_cells()
data = getter(self.cells) data = getter(self.cells)
if logbins: if logbins:
@@ -786,7 +788,7 @@ class SnapshotProcessor(HDF5Container):
if self.params.process.unload_cells: if self.params.process.unload_cells:
self.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]) centers = 0.5 * (edges[1:] + edges[:-1])
return (np.stack([values, centers]), {"logbins": logbins}) return (np.stack([values, centers]), {"logbins": logbins})
@@ -804,7 +806,7 @@ class SnapshotProcessor(HDF5Container):
self.unload_cells() self.unload_cells()
return ({"rho": centers, "B": B_mean}, {"logbins": logbins}) 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 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"]) dot_prod = np.einsum("ij,ij->i", vel_fluct, dset["Br"])
return np.abs(dot_prod) / (v_norm * B_norm) 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"]): def _mwa_sigma(self, axes=["x", "y", "z"]):
mw_speed = self.get_value("/globals/mwa_speed") mw_speed = self.get_value("/globals/mwa_speed")
@@ -905,7 +907,7 @@ class SnapshotProcessor(HDF5Container):
sigma_squared = sigma_squared + sigma_sq_ax sigma_squared = sigma_squared + sigma_sq_ax
return sigma_squared 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): def _coldens(self, ax_los):
datamap = self._rt.process(self._cam[ax_los], surf_qty=True) 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) dmap = self.get_value("/maps/" + name + "_" + ax_los)
if mass_weighted: 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: 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 return dmap / avg_map
def _dispersion_rad(self, name, ax_los="z"): 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 # 1. Get full storage names for the map and its azimuthal avererage
map_name = f"/maps/{name}_{ax_los}" 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 # 2. Get maps from the storage
dmap = self.get_value(map_name) dmap = self.get_value(map_name)
@@ -1236,7 +1238,7 @@ class SnapshotProcessor(HDF5Container):
return np.stack([values, centers]) return np.stack([values, centers])
def _fit_pdf(self, name, ax_los="z"): 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 values, centers = pdf
mask_fit = ( mask_fit = (
(centers > self.params.pdf.xmin_fit) (centers > self.params.pdf.xmin_fit)
@@ -1260,7 +1262,7 @@ class SnapshotProcessor(HDF5Container):
# Mean part # 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) radial_bins = self.get_value("/radial/radial_bins_" + ax_los)
mean_bin_vr = self.get_value("/radial/rad_avg_" + "slice_velr" + "_" + ax_los) mean_bin_vr = self.get_value("/radial/rad_avg_" + "slice_velr" + "_" + ax_los)
@@ -1299,7 +1301,7 @@ class SnapshotProcessor(HDF5Container):
return alpha return alpha
alpha_f = ( 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 # alpha
@@ -1310,8 +1312,8 @@ class SnapshotProcessor(HDF5Container):
"Map of the gravitational contribution to the Shakura&Sunaev alpha parameter for disks" "Map of the gravitational contribution to the Shakura&Sunaev alpha parameter for disks"
assert ax_los == "z" assert ax_los == "z"
T_avg = self.get_value("/maps/avg_map_T_mwavg_z") T_avg = self.get_value("/maps/rad_avg_map_T_mwavg_z")
coldens = self.get_value("/maps/avg_map_coldens_z") coldens = self.get_value("/maps/rad_avg_map_coldens_z")
def getter_alpha_grav(dset): def getter_alpha_grav(dset):
r2 = np.sum((self.lbox * self.oct_getter_pos_disk(dset)) ** 2, axis=2) 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") gphi = self.oct_getter_vect_phi(dset, "g")
return gr * gphi / (4 * np.pi * self.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 coldens * T_avg
) )
@@ -1515,6 +1517,70 @@ class SnapshotProcessor(HDF5Container):
return {"gfil": gfil, "Rfil": Rfil, "fPfil": fPfil, "dvr": dvr_fil} 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): def def_rules(self):
self.rules = { self.rules = {
@@ -1529,7 +1595,7 @@ class SnapshotProcessor(HDF5Container):
"T_mwavg": Rule( "T_mwavg": Rule(
self, self,
partial( partial(
self._ax_avg, self.ax_avg,
getter_T, getter_T,
mass_weighted=True, mass_weighted=True,
unit=self.info["unit_pressure"] / self.info["unit_density"], unit=self.info["unit_pressure"] / self.info["unit_density"],
@@ -1702,7 +1768,7 @@ class SnapshotProcessor(HDF5Container):
# PDF # PDF
"rho_pdf": Rule( "rho_pdf": Rule(
self, self,
partial(self._vol_pdf, partial(simple_getter, "rho"), logbins=True), partial(self.vol_pdf, partial(simple_getter, "rho"), logbins=True),
"Global rho-PDF", "Global rho-PDF",
"/hist", "/hist",
unit=self.info["unit_density"], unit=self.info["unit_density"],
@@ -1710,7 +1776,7 @@ class SnapshotProcessor(HDF5Container):
"rho_pdf_mw": Rule( "rho_pdf_mw": Rule(
self, self,
partial( partial(
self._vol_pdf, self.vol_pdf,
partial(simple_getter, "rho"), partial(simple_getter, "rho"),
weight_func=mass_func, weight_func=mass_func,
logbins=True, logbins=True,
@@ -1721,18 +1787,11 @@ class SnapshotProcessor(HDF5Container):
), ),
"T_pdf": Rule( "T_pdf": Rule(
self, self,
partial(self._vol_pdf, getter_T, logbins=True), partial(self.vol_pdf, getter_T, logbins=True),
"Global T-PDF", "Global T-PDF",
"/hist", "/hist",
unit=self.info["unit_pressure"] / self.info["unit_density"], 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( "cos_pdf": Rule(
self, self,
partial(self.cos_vfluct_B), partial(self.cos_vfluct_B),
@@ -1782,14 +1841,14 @@ class SnapshotProcessor(HDF5Container):
), ),
"mass": Rule( "mass": Rule(
self, self,
partial(self._sum, mass_func, mass_weighted=False), partial(self.sum, mass_func, mass_weighted=False),
"Total mass", "Total mass",
"/globals", "/globals",
unit=self.info["unit_density"] * self.info["unit_length"] ** 3, unit=self.info["unit_density"] * self.info["unit_length"] ** 3,
), ),
"mwa_speed": Rule( "mwa_speed": Rule(
self, self,
partial(self._vol_avg, partial(simple_getter, "vel")), partial(self.vol_avg, partial(simple_getter, "vel")),
"Mass weighted speed average", "Mass weighted speed average",
"/globals", "/globals",
unit=self.info["unit_velocity"], unit=self.info["unit_velocity"],
@@ -1804,14 +1863,14 @@ class SnapshotProcessor(HDF5Container):
), ),
"mwa_B_int": Rule( "mwa_B_int": Rule(
self, self,
partial(self._vol_avg, getter_B_int), partial(self.vol_avg, getter_B_int),
"Mass weighted Magnetic intensity average", "Mass weighted Magnetic intensity average",
"/globals", "/globals",
unit=self.info["unit_mag"], unit=self.info["unit_mag"],
), ),
} }
averageables = [ averageable_maps = [
"coldens", "coldens",
"Q", "Q",
"T", "T",
@@ -1829,7 +1888,7 @@ class SnapshotProcessor(HDF5Container):
self.rules["slice_" + name] = Rule( self.rules["slice_" + name] = Rule(
self, self,
partial(self._slice, getter, z=0.0, unit=unit), partial(self.slice, getter, z=0.0, unit=unit),
"{} slice".format(name), "{} slice".format(name),
"/maps", "/maps",
unit=unit, unit=unit,
@@ -1837,7 +1896,7 @@ class SnapshotProcessor(HDF5Container):
self.rules[name + "_mwavg"] = Rule( self.rules[name + "_mwavg"] = Rule(
self, 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), "Ax mass-weighted averaged {}".format(name),
"/maps", "/maps",
unit=unit, unit=unit,
@@ -1845,15 +1904,15 @@ class SnapshotProcessor(HDF5Container):
self.rules[name + "_avg"] = Rule( self.rules[name + "_avg"] = Rule(
self, 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), "Ax averaged {}".format(name),
"/maps", "/maps",
unit=unit, unit=unit,
) )
averageables.append("slice_" + name) averageable_maps.append("slice_" + name)
averageables.append(name + "_mwavg") averageable_maps.append(name + "_mwavg")
averageables.append(name + "_avg") averageable_maps.append(name + "_avg")
# special for vectors # special for vectors
if field in ["g", "vel", "Bl", "Br"]: 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]) generic_rule(field, partial(simple_getter, field), self.unit_key[field])
# radial average and other # radial average and other
for name in averageables: for name in averageable_maps:
unit = self.rules[name].unit unit = self.rules[name].unit
self.rules["rad_avg_" + name] = Rule( self.rules["rad_avg_" + name] = Rule(
self, self,
partial(self._rad_avg, name), partial(self._rad_avg, name),
@@ -1911,7 +1972,7 @@ class SnapshotProcessor(HDF5Container):
dependencies=["coldens", "radial_centers", "bins_on_map", name], dependencies=["coldens", "radial_centers", "bins_on_map", name],
unit=unit, unit=unit,
) )
self.rules["avg_map_" + name] = Rule( self.rules["rad_avg_map_" + name] = Rule(
self, self,
partial(self._rad_avg_map, name), partial(self._rad_avg_map, name),
"Interpolated map of azimuthal average of {}".format(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], dependencies=["radial_centers", "bins_on_map", "rad_avg_" + name],
unit=unit, unit=unit,
) )
self.rules["mwavg_map_" + name] = Rule( self.rules["rad_mwavg_map_" + name] = Rule(
self, self,
partial(self._rad_avg_map, name, mass_weighted=True), partial(self._rad_avg_map, name, mass_weighted=True),
"Interpolated map of azimuthal average of {}".format(name), "Interpolated map of azimuthal average of {}".format(name),
@@ -1932,14 +1993,14 @@ class SnapshotProcessor(HDF5Container):
partial(self._fluct_map, name), partial(self._fluct_map, name),
"Fluctuation wrt to average of {}".format(name), "Fluctuation wrt to average of {}".format(name),
"/maps", "/maps",
dependencies=[name, "avg_map_" + name], dependencies=[name, "rad_avg_map_" + name],
) )
self.rules["dispersion_rad_" + name] = Rule( self.rules["dispersion_rad_" + name] = Rule(
self, self,
partial(self._dispersion_rad, name), partial(self._dispersion_rad, name),
"radial RMS of {}".format(name), "radial RMS of {}".format(name),
"/radial", "/radial",
dependencies=[name, "avg_map_" + name], dependencies=[name, "rad_avg_map_" + name],
unit=unit, unit=unit,
) )
self.rules["mwfluct_" + name] = Rule( self.rules["mwfluct_" + name] = Rule(
@@ -1947,16 +2008,16 @@ class SnapshotProcessor(HDF5Container):
partial(self._fluct_map, name, mass_weigthed=True), partial(self._fluct_map, name, mass_weigthed=True),
"Fluctuation wrt to mass weighted average of {}".format(name), "Fluctuation wrt to mass weighted average of {}".format(name),
"/maps", "/maps",
dependencies=[name, "mwavg_map_" + name], dependencies=[name, "rad_mwavg_map_" + name],
) )
self.rules["pdf_" + name] = Rule( self.rules["fluct_pdf_" + name] = Rule(
self, self,
partial(self._rad_fluct_pdf, name), partial(self._rad_fluct_pdf, name),
"Probability density function of {} fluctuations".format(name), "Probability density function of {} fluctuations".format(name),
"/hist", "/hist",
dependencies=["rr", "fluct_" + name], dependencies=["rr", "fluct_" + name],
) )
self.rules["mwpdf_" + name] = Rule( self.rules["fluct_mwpdf_" + name] = Rule(
self, self,
partial(self._rad_fluct_pdf, name, mass_weigthed=True), partial(self._rad_fluct_pdf, name, mass_weigthed=True),
"Probability density function of {} mass weighted fluctuations".format( "Probability density function of {} mass weighted fluctuations".format(
@@ -1965,19 +2026,19 @@ class SnapshotProcessor(HDF5Container):
"/hist", "/hist",
dependencies=["rr", "mwfluct_" + name], dependencies=["rr", "mwfluct_" + name],
) )
self.rules["fit_pdf_" + name] = Rule( self.rules["fit_fluct_pdf_" + name] = Rule(
self, self,
partial(self._fit_pdf, name), partial(self._fit_pdf, name),
"Fit the PDF of {} fluctuations".format(name), "Fit the PDF of {} fluctuations".format(name),
"/hist", "/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 unit_bin = self.rules[name_bin].unit
if name_bin is not name: if name_bin is not name:
self.rules["mbb_" + name + "_" + name_bin] = Rule( self.rules["mbb_" + name + "_" + name_bin] = Rule(
self, 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), "Mean of {} by bins of {}".format(name, name_bin),
"/hist", "/hist",
dependencies=[name, name_bin], dependencies=[name, name_bin],