# coding: utf-8 import numpy as np import tables import pickle import astropy.units as u import pandas as pd from skimage.morphology import medial_axis import os from functools import partial from scipy.stats import linregress from astrophysix.simdm.results import Snapshot import pymses import pymses.utils.regions as reg from pymses.analysis import ( Camera, FractionOperator, MaxLevelOperator, ScalarOperator, raytracing, slicing, splatting, ) from pymses.filters import CellsToPoints, RegionFilter from fil_finder import FilFinder2D import pspec_new from units import U from baseprocessor import HDF5Container, Rule, norm_getter, simple_getter, vect_getter # Getters def mass_func(dset): dx = dset["dx"] return dset["rho"] * dx ** 3 # Mass function def vol_func(dset): return dset["dx"] ** 3 # Volume function def getter_T(dset): return dset["P"] / dset["rho"] # Temperature def getter_P(dset): return dset["P"] def getter_abs_cos_vB(dset): B_norm = np.sqrt(np.sum(dset["Br"] ** 2, axis=1)) v_norm = np.sqrt(np.sum(dset["vel"] ** 2, axis=1)) # Compute the dot product in each cell dot_prod = np.einsum("ij,ij->i", dset["vel"], dset["Br"]) return np.abs(dot_prod) / (v_norm * B_norm) def getter_B_int(dset): B_norm = np.sqrt(np.sum(dset["Br"] ** 2, axis=1)) return B_norm def getter_rho(dset): return dset["rho"] def getter_v_norm(dset): v_norm = np.sqrt(np.sum(dset["vel"] ** 2, axis=1)) return v_norm # Helpers def mean_by_bins( x, y, bins=100, logbins=False, ): """ Compute the mean of y in bins of x Parameters ---------- x, y : np.array of same dimensions bins : int, number of bins logbins : bool, if true, the bins will be logaritmically distributed """ mask = np.isfinite(x) & np.isfinite(y) x = x[mask].flatten() y = y[mask].flatten() if logbins: minvalue = np.min(x[x > 0]) x_bins = np.logspace(np.log10(minvalue), np.log10(np.max(x)), bins, base=10) else: x_bins = np.linspace(np.min(x), np.max(x), bins) # For each cell, bin_number contains the number of the bins it belongs to bin_number = np.zeros(len(y)) # Go through the min value of x of each bin for x_min in x_bins[1:-1]: bin_number = bin_number + (x > x_min).astype(int) # Compute the mean in each bin y_mean = np.zeros(len(x_bins) - 1) for i in range(len(y_mean)): y_mean[i] = np.mean(y[bin_number == i]) # Get the center of each bin if logbins: centers = 10 ** (0.5 * (np.log10(x_bins[1:]) + np.log10(x_bins[:-1]))) else: centers = 0.5 * (x_bins[1:] + x_bins[:-1]) return centers, y_mean # Filament helpers def find_center(distance, skeleton, i_center, j_center, i, j): """ Given a distance array, find the cells at a center of a filament at a given postion """ if skeleton[i, j]: i_center[i, j], j_center[i, j] = i, j return i, j elif i_center[i, j] or j_center[i, j]: return i_center[i, j], j_center[i, j] else: i_neigh = np.array([i - 1, i, i + 1]) i_neigh = i_neigh[(i_neigh > 0) & (i_neigh < distance.shape[0])] j_neigh = np.array([j - 1, j, j + 1]) j_neigh = j_neigh[(j_neigh > 0) & (j_neigh < distance.shape[1])] ii_neigh, jj_neigh = np.meshgrid(i_neigh, j_neigh) d_neigh = distance[ii_neigh, jj_neigh] ind_max = np.unravel_index(np.argmax(d_neigh), d_neigh.shape) i_max, j_max = ii_neigh[ind_max], jj_neigh[ind_max] if i_max == i and j_max == j: i_center[i, j], j_center[i, j] = i, j else: i_center[i, j], j_center[i, j] = find_center( distance, skeleton, i_center, j_center, i_max, j_max ) return i_center[i, j], j_center[i, j] # PostProcessor class class PostProcessor(HDF5Container): """ This class enable to compute and save derived quantities from the raw output """ # Axes information _ax_nb = {"x": 0, "y": 1, "z": 2} # Number of each axes _axes_h = {"x": "y", "y": "x", "z": "x"} # Associated horizontal axe _axes_v = {"x": "z", "y": "z", "z": "y"} # Associated vertical axe # Pymses unit key of amr fiels unit_key = { "rho": "unit_density", "vel": "unit_velocity", "Br": "unit_mag", "Bl": "unit_mag", "P": "unit_pressure", "g": {"unit_gravpot": 1, "unit_length": -1}, "phi": "unit_gravpot", } G = 1.0 # Gravitational constant cells_loaded = False parts_loaded = False def __init__( self, path=None, num=None, path_out=None, pp_params=None, tag=None, unit_time=U.year, ): """ Creates the basic structures needed for the outputs """ super(PostProcessor, self).__init__(path, path_out, pp_params, tag) # Open outfile if not self.pp_params.out.tag == "": tag_name = self.pp_params.out.tag + "_" else: tag_name = "" self.filename = ( self.path_out + "/postproc_" + tag_name + format(num, "05") + ".h5" ) self.cells_filename = ( self.path_out + "/cells_" + tag_name + format(num, "05") + ".h5" ) self.parts_filename = ( self.path_out + "/parts_" + tag_name + format(num, "05") + ".h5" ) self.pspec_filename = ( self.path_out + "/pspec_" + tag_name + format(num, "05") + ".h5" ) self.filaments_filename = ( self.path_out + "/filaments_" + tag_name + format(num, "05") + ".pickle" ) if not os.path.exists(self.path_out): os.makedirs(self.path_out) self.open() # Ramses Output self.path = path self.run = os.path.basename(path) self.num = num self._ro = pymses.RamsesOutput( path, num, order=self.pp_params.pymses.order, verbose=self.pp_params.pymses.verbose, ) self._amr = self._ro.amr_source(self.pp_params.pymses.variables) self._part = self._ro.particle_source(["vel", "mass"]) if self.pp_params.pymses.filter: self.min_coords = np.array(self.pp_params.pymses.min_coords) self.max_coords = np.array(self.pp_params.pymses.max_coords) box = reg.Box((self.min_coords, self.max_coords)) amr_filt = RegionFilter(box, self._amr) self._amr = amr_filt part_filt = RegionFilter(box, self._part) self._part = part_filt self.info = self._ro.info.copy() # Density operator self._rho_op = ScalarOperator( lambda dset: dset["rho"], self._ro.info["unit_density"] ) # Density ray tracer if self.pp_params.pymses.fft: self._rt = splatting.SplatterProcessor( self._amr, self._ro.info, self._rho_op ) else: self._rt = raytracing.RayTracer(self._amr, self._ro.info, self._rho_op) if not self.pp_params.pymses.multiprocessing: self._rt.disable_multiprocessing() # Set the extent of the image self._radius = 0.5 / self.pp_params.pymses.zoom self.lbox = self.info["boxlen"] if self.pp_params.pymses.filter: center = (self.max_coords + self.min_coords) / 2.0 im_extent = np.array( [ self.min_coords[0], self.max_coords[0], self.min_coords[1], self.max_coords[1], ] ) distance = (self.max_coords[2] - self.min_coords[2]) / 2.0 else: center = self.pp_params.pymses.center im_extent = np.array( [ (-self._radius + center[0]), (self._radius + center[0]), (-self._radius + center[1]), (self._radius + center[1]), ] ) distance = self._radius # Get time time = self._ro.info["time"] # time in codeunits # Set post processing attributes self.save.root._v_attrs.dir = os.path.dirname(path) self.save.root._v_attrs.run = os.path.basename(path) self.save.root._v_attrs.num = num self.save.root._v_attrs.lbox = self.lbox self.save.root._v_attrs.unit_length = self.info["unit_length"] self.save.root._v_attrs.time = time if "/maps" not in self.save: self.save.create_group("/", "maps", "2D maps") self.save.root.maps._v_attrs.center = center self.save.root.maps._v_attrs.radius = self._radius self.save.root.maps._v_attrs.im_extent = im_extent # Initialize cameras self._cam = {} for ax_los in self._ax_nb: # los = line of sight ax_v = self._axes_v[ax_los] self._cam[ax_los] = Camera( center=self.pp_params.pymses.center, line_of_sight_axis=ax_los, region_size=[im_extent[1] - im_extent[0], im_extent[3] - im_extent[2]], distance=distance, far_cut_depth=distance, up_vector=ax_v, map_max_size=self.pp_params.pymses.map_size, ) self.close() self.log_id = "[{}, {}] ".format(self.run, self.num) if os.path.exists(self.filaments_filename): with open(self.filaments_filename, "rb") as f: self.fil = pickle.load(f) else: self.fil = None time_in_right_unit = self.info["time"] * self.info["unit_time"].express( unit_time ) self.snapshot = Snapshot( name=str(self.num), description="", time=(time_in_right_unit, unit_time), directory_path=self.path, data_reference="OUTPUT_{}".format(self.num), ) self.def_rules() def load_data(self, points_src, loaded, filename, save): """ Load data from the source file in the memory. (Long and memory heavy) """ if not self.parts_loaded: if os.path.exists(filename): hdf5 = tables.open_file(filename, mode="r") try: node = hdf5.get_node("/data") data = {} for key in node._v_children: data[key] = hdf5.get_node("/data/" + key).read() finally: hdf5.close() else: data_pymses = points_src.flatten() data = {} for key in data_pymses.fields: data[key] = data_pymses[key] try: data["dx"] = data_pymses.get_sizes() except AttributeError: pass data["pos"] = data_pymses.points if save: hdf5 = tables.open_file(filename, mode="w") try: for key in data: hdf5.create_array( "/data", key, data[key], "", createparents=True ) finally: hdf5.close() return data def load_parts(self): self.parts = self.load_data( self._part, self.parts_loaded, self.parts_filename, self.pp_params.process.save_parts, ) self.parts_loaded = True def unload_parts(self): """ Free space in the memory by telling the garbage collectors that self.parts is not needed """ if self.parts_loaded: del self.parts self.parts_loaded = False def load_cells(self): """ Load all cells from the source file in the memory. Cells will be accessible trough self.cells (Long and memory heavy) """ cells_src = CellsToPoints(self._amr) self.cells = self.load_data( cells_src, self.cells_loaded, self.cells_filename, self.pp_params.process.save_cells, ) self.parts_loaded = True def unload_cells(self): """ Free space in the memory by telling the garbage collectors that self.cells is not needed """ if self.cells_loaded: del self.cells self.cells_loaded = False def getter_pos_disk(self, dset): """ Returns the position in normalized and centered units """ pos = dset.points pos = pos - np.array(self.save.root.maps._v_attrs.center) return pos def getter_vect_r(self, dset, name_vect): """ Radial component of a vector """ r = self.getter_pos_disk(dset)[:, :2] ur = np.transpose((np.transpose(r) / np.sqrt(np.sum(r ** 2, axis=1)))) return np.einsum("ij, ij -> i", dset[name_vect][:, :2], ur) def getter_vect_phi(self, dset, name_vect): """ Azimuthal component of a vector """ r = self.getter_pos_disk(dset)[:, :2] r_norm = np.sqrt(np.sum(r ** 2, axis=1)) rot = np.array([[0, -1], [1, 0]]) uphi = np.transpose(np.einsum("ij, kj -> ik", rot, r) / r_norm) vect = dset[name_vect][:, :2] return np.einsum("ij,ij -> i", vect, uphi) def oct_getter_pos_disk(self, dset): """ Returns the position in normalized units centered on the position of the star """ pos = dset.get_cell_centers() pos = pos - (np.array(self.pp_params.disk.pos_star) / self.lbox) return pos def oct_getter_vect_r(self, dset, name_vect): """ Radial component of a vector """ r = self.oct_getter_pos_disk(dset)[:, :, :2] ur = np.transpose( (np.transpose(r, (2, 0, 1)) / np.sqrt(np.sum(r ** 2, axis=2))), (1, 2, 0) ) return np.einsum("ikj, ikj -> ik", dset[name_vect][:, :, :2], ur) def oct_getter_vect_phi(self, dset, name_vect): """ Azimuthal component of a vector """ r = self.oct_getter_pos_disk(dset)[:, :, :2] r_norm = np.sqrt(np.sum(r ** 2, axis=2)) rot = np.array([[0, -1], [1, 0]]) uphi = np.transpose(np.einsum("ij, klj -> ikl", rot, r) / r_norm, (1, 2, 0)) vect = dset[name_vect][:, :, :2] return np.einsum("ikj,ikj -> ik", vect, uphi) def oct_getter_vr(self, dset): return self.oct_getter_vect_r(dset, "vel") def oct_getter_vphi(self, dset): """ Azimuthal velocity """ return self.oct_getter_vect_phi(dset, "vel") def _slice(self, getter, ax_los="z", z=0.0, unit=U.none): """ Slice process function. Return a slice of the source box. Parameters ---------- getter : callable A callable that extract the wanted data from a pymses dataset ax_los : string The axis perpendicular to the slice plane z : float Coordinate of the slice on the ax_los axis unit : U.Unit Unit of the resulting dataset Returns ------- A numpy array containing the slice """ unit = self._get_units(unit) op = ScalarOperator(getter, unit) 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): """ Map of the average of a quantity (given by getter) along an axis (ax_los) Returns 2D array if getter returns a scalar quantity If surf_qty is set (projection mode), mass_weighted is ignored """ unit = self._get_units(unit) if surf_qty: op = ScalarOperator(getter, unit) else: if mass_weighted: def num(dset): value = getter(dset) rho = getter_rho(dset) return rho * value def denum(dset): return getter_rho(dset) else: # Volume weighted def num(dset): value = getter(dset) return value def denum(dset): return 1.0 op = FractionOperator(num, denum, unit) if self.pp_params.pymses.fft: rt = splatting.SplatterProcessor(self._amr, self._ro.info, op) else: rt = raytracing.RayTracer(self._amr, self._ro.info, op) if not self.pp_params.pymses.multiprocessing: rt.disable_multiprocessing() datamap = rt.process(self._cam[ax_los], surf_qty=surf_qty) return datamap.map.T def _get_axis(self, axis): if isinstance(axis, str): axis = self._ax_nb[axis] self.load_cells() return np.sort(np.unique(self.cells["pos"][:, axis])) def _plane_avg_uniform(self, getter, axis, unit=U.none, mass_weighted=False): """ Profile of the average of a quantity (given by getter) perpendicular to an axis WARNING : This version only works on an uniform grid, need of a box version for AMR Returns 1D array if getter returns a scalar quantity """ unit = self._get_units(unit) self.load_cells() if isinstance(axis, str): axis = self._ax_nb[axis] axis_data = self.cells["pos"][:, axis] value = getter(self.cells) df = pd.DataFrame({"axis": axis_data}) if mass_weighted: mass = mass_func(self.cells) tot_mass = np.sum(mass) df["value"] = value * mass / tot_mass else: df["value"] = value if self.pp_params.process.unload_cells: self.unload_cells() df.sort_values("axis", inplace=True) return df.groupby("axis").mean().values[:, 0] 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) """ self.load_cells() value = getter(self.cells) if mass_weighted: mass = mass_func(self.cells) # Transpose (.T) is for vectorial values data = np.sum((mass * value.T).T, axis=0) / np.sum(mass) else: data = np.sum(value, axis=0) if self.pp_params.process.unload_cells: self.unload_cells() return data 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) """ self.load_cells() value = getter(self.cells) if mass_weighted: weight = mass_func(self.cells) else: weight = vol_func(self.cells) # Transpose (.T) is for vectorial values data = np.sum((weight * value.T).T, axis=0) / np.sum(weight) if self.pp_params.process.unload_cells: self.unload_cells() return data def _vol_pdf(self, getter, bins=100, logbins=False, weight_func=vol_func): self.load_cells() data = getter(self.cells) if logbins: data = np.log10(data) weights = weight_func(self.cells) if self.pp_params.process.unload_cells: self.unload_cells() values, edges = np.histogram(data, bins, weights=weights) centers = 0.5 * (edges[1:] + edges[:-1]) return (np.stack([values, centers]), {"logbins": logbins}) def _Brho(self, bins=100, logbins=True): """ Mean of B in rho bins """ self.load_cells() B = getter_B_int(self.cells) rho = getter_rho(self.cells) centers, B_mean = mean_by_bins(rho, B, bins, logbins) if self.pp_params.process.unload_cells: self.unload_cells() return ({"rho": centers, "B": B_mean}, {"logbins": logbins}) def _mean_by_bins( self, name_x, name_y, ax_los, group="/maps/", bins=100, logbins=True ): """ Compute the mean of y in bins of x, where x, y are to array of the hdf5 file Parameters ---------- x_name, y_name : str, path of x and y in the hdf5 file bins : int, number of bins logbins : bool, if true, the bins will be logaritmically distributed """ x = self.save.get_node(group + name_x + "_" + ax_los).read() y = self.save.get_node(group + name_y + "_" + ax_los).read() centers, values = mean_by_bins(x, y, bins, logbins) # return ({os.path.basename(name_x) : centers, os.path.basename(name_y) : y_mean}, # {"logbins" : logbins}) return (np.stack([values, centers]), {"logbins": logbins}) def _Ek_Eb_rho(self, bins=100, logbins=True): """ Mean of Ek/Eb in rho bins """ self.load_cells() mean_speed = self.save.get_node("/globals/mwa_speed").read() mean_speed = mean_speed * self.info["unit_velocity"].express(U.km_s) vel_fluct = (self.cells)["vel"] * self.info["unit_velocity"].express( U.km_s ) - mean_speed B_norm = getter_B_int(self.cells) B_norm = B_norm * self.info["unit_mag"].express(U.T) v_norm = np.sqrt( np.sum((vel_fluct * 10 ** (3)) ** 2, axis=1) ) # v_norm [m/s] et vel_fluct [km/s] rho = getter_rho(self.cells) rho_kg_m3 = rho * self.info["unit_density"].express(U.kg_m3) eb = 0.5 * (B_norm) ** 2 / (4 * np.pi * 10 ** (-7)) # mettre le bon mu ek = 0.5 * v_norm ** 2 * rho_kg_m3 rapport = ek / eb if logbins: rho_bins = np.logspace( np.log10(np.min(rho)), np.log10(np.max(rho)), bins, base=10 ) else: rho_bins = np.linspace(np.min(rho), np.max(rho), bins) # For each cell, bin_number contains the number of the bins it belongs to bin_number = np.zeros(len(B_norm)) # Go through the min value of rho of each bin for rho_min in rho_bins[:-1]: bin_number = bin_number + (rho > rho_min).astype(int) # Compute the mean in each bin ek_eb = np.zeros(len(rho_bins) - 1) for i in range(len(ek_eb)): ek_eb[i] = np.mean(rapport[bin_number == i]) # Get the center of each bin if logbins: centers = 10 ** (0.5 * (np.log10(rho_bins[1:]) + np.log10(rho_bins[:-1]))) else: centers = 0.5 * (rho_bins[1:] + rho_bins[:-1]) if self.pp_params.process.unload_cells: self.unload_cells() return ({"rho": centers, "Ek_Eb_rho": ek_eb}, {"logbins": logbins}) def cos_vfluct_B(self): "return the cos of the angle between the magnetic field and the velocity fluctuation field" mean_speed = self.save.get_node("/globals/mwa_speed").read() def getter_cos_vfluct_B(dset): vel_fluct = dset["vel"] - mean_speed B_norm = np.sqrt(np.sum(dset["Br"] ** 2, axis=1)) v_norm = np.sqrt(np.sum(vel_fluct ** 2, axis=1)) # Compute the dot product in each cell 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) def _mwa_sigma(self, axes=["x", "y", "z"]): mw_speed = self.save.get_node("/globals/mwa_speed").read() if axes == ["x", "y", "z"]: def getter(dset): return np.sum((dset["vel"] - mw_speed) ** 2, axis=1) else: def getter(dset): sigma_squared = 0.0 for ax in axes: ax_nb = self._ax_nb[ax] sigma_sq_ax = (dset["vel"][:, ax_nb] - mw_speed[ax_nb]) ** 2 sigma_squared = sigma_squared + sigma_sq_ax return sigma_squared 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) return datamap.map.T def _rho(self, ax_los, z=0.0): datamap_rho = slicing.SliceMap(self._amr, self._cam[ax_los], self._rho_op, z=z) return (datamap_rho.map).T def _vector_h(self, name, unit, ax_los, z=0.0): h_op = ScalarOperator( lambda dset: dset[name][:, self._ax_nb[self._axes_h[ax_los]]], self._ro.info[unit], ) dmap_h = slicing.SliceMap(self._amr, self._cam[ax_los], h_op, z=z).map.T return dmap_h def _vector_v(self, name, unit, ax_los, z=0.0): v_op = ScalarOperator( lambda dset: dset[name][:, self._ax_nb[self._axes_v[ax_los]]], self._ro.info[unit], ) dmap_v = slicing.SliceMap(self._amr, self._cam[ax_los], v_op, z=z).map.T return dmap_v def _speed_h(self, ax_los, z=0.0): return self._vector_h("vel", "unit_velocity", ax_los, z) def _speed_v(self, ax_los, z=0.0): return self._vector_v("vel", "unit_velocity", ax_los, z) def _B_h(self, ax_los, z=0.0): return self._vector_h("Br", "unit_mag", ax_los, z) def _B_v(self, ax_los, z=0.0): return self._vector_v("Br", "unit_mag", ax_los, z) def _B_int(self, ax_los, z=0.0): """ Slice ont the intensity of the magnetic field """ B_op = ScalarOperator( lambda dset: np.sqrt(np.sum(dset["Br"] ** 2, axis=1)), self._ro.info["unit_mag"], ) dmap_B = (slicing.SliceMap(self._amr, self._cam[ax_los], B_op, z=z)).map.T return dmap_B def _temperature(self, ax_los, z=0.0): P_op = ScalarOperator(lambda dset: dset["P"], self._ro.info["unit_pressure"]) dmap_P = (slicing.SliceMap(self._amr, self._cam[ax_los], P_op, z=z)).map.T dmap_rho = self.save.get_node("/maps/rho_{}".format(ax_los)).read() return dmap_P / dmap_rho def _levels(self, ax_los): self._amr.set_read_levelmax(self.pp_params.pymses.levelmax) level_op = MaxLevelOperator() rt_level = raytracing.RayTracer(self._amr, self._ro.info, level_op) if not self.pp_params.pymses.multiprocessing: rt_level.disable_multiprocessing() datamap = rt_level.process(self._cam[ax_los], surf_qty=True) return datamap.map.T def _jeans(self, ax_los): dmap_T = self.save.get_node("/maps/T_" + ax_los).read() dmap_rho = self.save.get_node("/maps/rho_" + ax_los).read() dmap_jeans = np.sqrt(np.pi * dmap_T / dmap_rho) return dmap_jeans def _jeans_ratio(self, ax_los): dmap_jeans = self.save.get_node("/maps/jeans_" + ax_los).read() dmap_levels = self.save.get_node("/maps/levels_" + ax_los).read() dmap_jeans_ratio = dmap_jeans * 2 ** (dmap_levels) return dmap_jeans_ratio def _toomreQ_disk(self, ax_los): """ Compute the Toomre Q parameter in a Keplerian disk """ # Operator to compute the angular speed times rho def omega_rho_func(dset): pos = self.oct_getter_pos_disk(dset) xx = pos[:, :, 0] yy = pos[:, :, 1] rc = np.sqrt(xx ** 2 + yy ** 2) # cylindrical radius vx = dset["vel"][:, :, 0] vy = dset["vel"][:, :, 1] omega_rho = 1.0 / rc ** 2 omega_rho = omega_rho * dset["rho"] vyx = vy * xx vxy = vx * yy omega_rho = omega_rho * (vyx - vxy) return omega_rho # Operator to compute the angular speed omega_op = FractionOperator( omega_rho_func, lambda dset: dset["rho"], 1.0 / self._ro.info["unit_time"] ) # Operator to compute the sound speed cs_op = FractionOperator( lambda dset: dset["P"], lambda dset: dset["rho"], self._ro.info["unit_velocity"], ) # Ray tracer for the angular speed rt_omega = raytracing.RayTracer(self._amr, self._ro.info, omega_op) # Ray tracer for the sound speed if self.pp_params.pymses.fft: rt_cs = splatting.SplatterProcessor( self._amr, self._ro.info, cs_op, surf_qty=False ) else: rt_cs = raytracing.RayTracer(self._amr, self._ro.info, cs_op) if not self.pp_params.pymses.multiprocessing: rt_omega.disable_multiprocessing() rt_cs.disable_multiprocessing() dmap_omega = rt_omega.process(self._cam[ax_los]) dmap_cs = rt_cs.process(self._cam[ax_los]) dmap_col = self.save.root.maps.coldens_z.read() map_Q = ( (self.lbox * dmap_cs.map.T) * dmap_omega.map.T / (np.pi * self.G * dmap_col) ) return map_Q def _radial_bins(self, ax_los="z"): """ Computes radial bins (for disk) """ pos_star = self.pp_params.disk.pos_star im_extent = np.array(self.save.root.maps._v_attrs.im_extent) * self.lbox # radius of the corner of the box plus a margin rad_of_box = ( np.sqrt( (im_extent[1] - pos_star[0]) ** 2 + (im_extent[3] - pos_star[1]) ** 2 ) + 0.1 ) bin_in = self.pp_params.disk.bin_in bin_out = self.pp_params.disk.bin_out nb_bin = self.pp_params.disk.nb_bin # radial bins if self.pp_params.disk.binning == "log": lrad_in = np.log10(bin_in) lrad_ext = np.log10(bin_out) rad_bins = np.logspace(lrad_in, lrad_ext, num=nb_bin) elif self.pp_params.disk.binning == "lin": rad_bins = np.linspace(bin_in, bin_out, num=nb_bin) # Add boundaries rad_bins = np.concatenate(([0.0], rad_bins, [rad_of_box])) return rad_bins def _radial_centers(self, ax_los="z"): radial_bins = self.save.get_node("/radial/radial_bins_" + ax_los).read() bin_centers = 0.5 * (radial_bins[1:] + radial_bins[:-1]) return bin_centers def _rr(self, ax_los="z"): """ Computes the radius from the center """ im_extent = np.array(self.save.root.maps._v_attrs.im_extent) * self.lbox map_size = self.pp_params.pymses.map_size pos_star = self.pp_params.disk.pos_star # Physical size of cells dx = (im_extent[1] - im_extent[0]) / map_size dy = (im_extent[3] - im_extent[2]) / map_size # Physical coordinates of the center of the cells x = np.linspace(im_extent[0], im_extent[1], map_size) + 0.5 * dx y = np.linspace(im_extent[2], im_extent[3], map_size) + 0.5 * dy xx, yy = np.meshgrid(x, y) # Physical radius rr = np.sqrt((xx - pos_star[0]) ** 2 + (yy - pos_star[1]) ** 2) return rr def _bins_on_map(self, ax_los="z"): rad_bins = self.save.get_node("/radial/radial_bins_" + ax_los).read() rr = self.save.get_node("/maps/rr_" + ax_los).read() # Find appropriate bin for each coordinate set bins = np.zeros(rr.shape, dtype=int) for r in rad_bins[1:]: bins = bins + (rr >= r).astype(int) return bins def _rad_avg(self, name, ax_los="z"): radial_bins = self.save.get_node("/radial/radial_bins_" + ax_los).read() bins_on_map = self.save.get_node("/maps/bins_on_map_" + ax_los).read() dmap = self.save.get_node("/maps/" + name + "_" + ax_los).read() # mean of all the cells in the bin mean_bin = np.zeros(len(radial_bins) - 1) for j in range(len(radial_bins) - 1): mean_bin[j] = np.mean(dmap[bins_on_map == j]) return mean_bin def _rad_avg_map(self, name, ax_los="z"): radial_bins = self.save.get_node("/radial/radial_bins_" + ax_los).read() bins_on_map = self.save.get_node("/maps/bins_on_map_" + ax_los).read() rr = self.save.get_node("/maps/rr_" + ax_los).read() mean_bin = self.save.get_node("/radial/rad_avg_" + name + "_" + ax_los).read() # Add value for border mean_bin = np.concatenate((mean_bin, [mean_bin[-1]])) rr_flat = rr.flatten() bins_on_map_flat = bins_on_map.flatten() # Compute the map azimuthally averaged # use linear interpolation to improve accuracy avg_flat = (radial_bins[bins_on_map_flat + 1] - rr_flat) * mean_bin[ bins_on_map_flat ] avg_flat = ( avg_flat + (rr_flat - radial_bins[bins_on_map_flat]) * mean_bin[bins_on_map_flat + 1] ) avg_flat = avg_flat / ( radial_bins[bins_on_map_flat + 1] - radial_bins[bins_on_map_flat] ) avg_map = np.reshape(avg_flat, rr.shape) return avg_map def _fluct_map(self, name, ax_los="z"): dmap = self.save.get_node("/maps/" + name + "_" + ax_los).read() avg_map = self.save.get_node("/maps/avg_map_" + name + "_" + ax_los).read() return dmap / avg_map def _rad_fluct_pdf(self, name, ax_los="z"): fluct_map = self.save.get_node("/maps/fluct_" + name + "_" + ax_los).read() rr = self.save.get_node("/maps/rr_" + ax_los).read() mask_pdf = ( (rr > self.pp_params.disk.rmin_pdf) & (rr < self.pp_params.disk.rmax_pdf) & (fluct_map > 0) ) values, edges = np.histogram( np.log10(fluct_map[mask_pdf].flatten()), self.pp_params.pdf.nb_bin, range=self.pp_params.pdf.range, density=True, ) centers = 0.5 * (edges[1:] + edges[:-1]) return np.stack([values, centers]) def _fit_pdf(self, name, ax_los="z"): pdf = self.save.get_node("/hist/pdf_" + name + "_" + ax_los) values, centers = pdf.read() mask_fit = ( (centers > self.pp_params.pdf.xmin_fit) & (centers < self.pp_params.pdf.xmax_fit) & (values > self.pp_params.pdf.fit_cut * np.max(values)) ) (slope, origin, correlation, _, stderr) = linregress( centers[mask_fit], np.log10(values[mask_fit]) ) pdf.attrs.slope = slope pdf.attrs.origin = origin pdf.attrs.correlation = correlation pdf.attrs.stderr = stderr pdf.attrs.var = np.var return True def _alpha_disk(self, ax_los="z"): "Map of the Rayleigh contribution to the Shakura&Sunaev alpha parameter for disks" assert ax_los == "z" # Mean part T_avg = self.save.get_node("/maps/avg_map_T_mwavg_z").read() radial_bins = self.save.get_node("/radial/radial_bins_" + ax_los).read() mean_bin_vr = self.save.get_node( "/radial/rad_avg_" + "vr" + "_" + ax_los ).read() mean_bin_vphi = self.save.get_node( "/radial/rad_avg_" + "vphi" + "_" + ax_los ).read() mean_bin_vr = np.concatenate((mean_bin_vr, [mean_bin_vr[-1]])) mean_bin_vphi = np.concatenate((mean_bin_vphi, [mean_bin_vr[-1]])) # Fluct part def getter_alpha_num(dset): r = np.sqrt( np.sum((self.lbox * self.oct_getter_pos_disk(dset)) ** 2, axis=2) ) bins = np.zeros(r.shape, dtype=int) for r0 in radial_bins[1:-1]: bins = bins + (r >= r0).astype(int) vr_mean = mean_bin_vr[bins] vphi_mean = mean_bin_vphi[bins] # use linear interpolation # v = ((r[i+1] - r)v[i] + (r - r[i])v[i + 1]) / (r[i+1] - r[i]) vr_mean = (radial_bins[bins + 1] - r) * vr_mean vr_mean = vr_mean + (r - radial_bins[bins]) * mean_bin_vr[bins + 1] vr_mean = vr_mean / (radial_bins[bins + 1] - radial_bins[bins]) vphi_mean = (radial_bins[bins + 1] - r) * vphi_mean vphi_mean = vphi_mean + (r - radial_bins[bins]) * mean_bin_vphi[bins + 1] vphi_mean = vphi_mean / (radial_bins[bins + 1] - radial_bins[bins]) vr = self.oct_getter_vr(dset) vphi = self.oct_getter_vphi(dset) alpha = (vphi - vphi_mean) * (vr - vr_mean) return alpha alpha_f = ( self._ax_avg(getter_alpha_num, "z", unit=U.none, mass_weighted=True) / T_avg ) # alpha alpha = (2.0 / 3) * alpha_f return alpha def _alpha_grav(self, ax_los="z"): "Map of the gravitational contribution to the Shakura&Sunaev alpha parameter for disks" assert ax_los == "z" T_avg = self.save.get_node("/maps/avg_map_T_mwavg_z").read() coldens = self.save.get_node("/maps/avg_map_coldens_z").read() def getter_alpha_grav(dset): r2 = np.sum((self.lbox * self.oct_getter_pos_disk(dset)) ** 2, axis=2) e2 = (1.0 / 256.0) ** 2 gstar = -self.G * self.pp_params.disk.mass_star / (e2 + r2) gr = self.oct_getter_vect_r(dset, "g") - gstar 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) / ( coldens * T_avg ) # alpha alpha_g = (2.0 / 3) * alpha_g return alpha_g def _sinks(self): csv_name = ( self.path + "/output_" + str(self.num).zfill(5) + "/sink_" + str(self.num).zfill(5) + ".csv" ) sinks = np.loadtxt(csv_name, delimiter=",") header = [ "Id", "M", "dmf", "x", "y", "z", "vx", "vy", "vz", "rot_period", "lx", "ly", "lz", "acc_rate", "acc_lum", "age", "int_lum", "Teff", ] if len(sinks) == 0: sinks = np.zeros(len(header)) sinks_dict = {} for key, a in zip(header, sinks): sinks_dict[key] = a return sinks_dict def _pspec(self): outfile = self.pspec_filename pspec_new.pspec(repo=self.path, iouts=[self.num], outfile=outfile) return True def _particles(self): """Ensure particles are written in the hdf5 file""" self.load_parts() self.unload_parts() return np.array([self.parts_filename]) def _filaments(self): datamap_name = self.pp_params.filaments.datamap verbose = self.pp_params.filaments.verbose rmin_frac = self.pp_params.filaments.rmin rmax_frac = self.pp_params.filaments.rmax size_thresh = self.pp_params.filaments.size_thresh skel_thresh = self.pp_params.filaments.skel_thresh branch_thresh = self.pp_params.filaments.branch_thresh glob_thresh = self.pp_params.filaments.glob_thresh datamap = self.save.get_node("/maps/" + datamap_name + "_z").read() shape = datamap.shape x = np.arange(shape[0]) - shape[0] / 2 y = np.arange(shape[1]) - shape[1] / 2 xx, yy = np.meshgrid(x, y) rr = np.sqrt(xx ** 2 + yy ** 2) rmin = int(rmin_frac * shape[0]) rmax = int(rmax_frac * shape[0]) mask = (rr >= rmin) & (rr <= rmax) datamap[np.logical_not(mask)] = np.nan self.fil = FilFinder2D(datamap, distance=1 * u.cm, beamwidth=1 * u.pix) self.fil.preprocess_image(flatten_percent=95) self.fil.create_mask( verbose=verbose, smooth_size=1 * u.pix, adapt_thresh=4 * u.pix, size_thresh=size_thresh * u.pix ** 2, glob_thresh=glob_thresh, fill_hole_size=0.1 * u.pix ** 2, ) self.fil.medskel(verbose=verbose) self.fil.analyze_skeletons( skel_thresh=skel_thresh * u.pix, branch_thresh=branch_thresh * u.pix, relintens_thresh=0.1, ) self.fil.exec_rht() self.fil.find_widths() with open(self.filaments_filename, "wb") as f: pickle.dump(self.fil, f, pickle.HIGHEST_PROTOCOL) return True def _filaments_center(self): """ Fill an array with center postion for each cell in a filament """ fil = self.fil mask = fil.mask.copy() _, distance = medial_axis(mask, return_distance=True) skel = fil.skeleton i_center = np.zeros(distance.shape, dtype=int) j_center = np.zeros(distance.shape, dtype=int) x_mask, y_mask = np.where(mask) for k in range(len(x_mask)): find_center(distance, skel, i_center, j_center, x_mask[k], y_mask[k]) return np.stack([i_center, j_center]) def _filaments_forces(self): """ Compute forces within a filament (for disks), within the slice at z=0 """ GM = self.G * self.pp_params.disk.mass_star # Mass parameter # Find center of filaments i_center, j_center = self._filaments_center() # Get slices and projections at z = 0 vphi = self.get_value("/maps/slice_velphi_z") gr = self.get_value("/maps/slice_gr_z") Pz = self.get_value("/maps/slice_P_z") rho = self.get_value("/maps/slice_rho_z") vr = self.get_value("/maps/slice_velr_z") # Get coordinates im_extent = np.array(self.save.root.maps._v_attrs.im_extent) * self.lbox map_size = self.pp_params.pymses.map_size pos_star = self.pp_params.disk.pos_star # Physical size of cells dx = (im_extent[1] - im_extent[0]) / map_size dy = (im_extent[3] - im_extent[2]) / map_size # Physical coordinates of the center of the cells x = np.linspace(im_extent[0], im_extent[1], map_size) + 0.5 * dx y = np.linspace(im_extent[2], im_extent[3], map_size) + 0.5 * dy xx, yy = np.meshgrid(x, y) rr = np.sqrt((xx - pos_star[0]) ** 2 + (yy - pos_star[1]) ** 2) # Rotational support R = vphi ** 2 / rr # Equilibrium Gvrx, Gvry = np.gradient(vr) gradvr = (xx * Gvrx + yy * Gvry) / rr dvr = gradvr + vr * gradvr # Complete derivative # Thermal support GPx, GPy = np.gradient(Pz) gradPr = (xx * GPx + yy * GPy) / rr fP = gradPr / rho # Gravitational field e2 = (1.0 / 512) ** 2 gstar = -GM / (rr ** 2 + e2) # Substract gravitational field from the star Rdisk = R + gstar gdisk = gr - gstar # Forces at the center of filaments Rdisk_center = Rdisk[i_center, j_center] gr_center = gdisk[i_center, j_center] fP_center = fP[i_center, j_center] dvr_center = dvr[i_center, j_center] # Forces for the filaments equilibrium Rfil = Rdisk - Rdisk_center gfil = gdisk - gr_center fPfil = fP - fP_center dvr_fil = dvr - dvr_center return {"gfil": gfil, "Rfil": Rfil, "fPfil": fPfil, "dvr": dvr_fil} def def_rules(self): self.rules = { # Base rules "coldens": Rule( self, self._coldens, "Column density", "/maps", unit=self.info["unit_density"] * self.info["unit_length"], ), "vr": Rule( self, partial( self._ax_avg, self.oct_getter_vr, mass_weighted=True, unit=self.info["unit_velocity"], ), "Vertically mass-weighted averaged radial velocity", "/maps", unit=self.info["unit_velocity"], ), "vphi": Rule( self, partial( self._ax_avg, self.oct_getter_vphi, mass_weighted=True, unit=self.info["unit_velocity"], ), "Vertically mass-weighted averaged azimuthal velocity", "/maps", unit=self.info["unit_velocity"], ), "T_mwavg": Rule( self, partial( self._ax_avg, getter_T, mass_weighted=True, unit=self.info["unit_pressure"] / self.info["unit_density"], ), "Ax mass-weighted averaged azimuthal temperature", "/maps", unit=self.info["unit_pressure"] / self.info["unit_density"], ), "alpha_disk": Rule( self, self._alpha_disk, "Map of the Shakura&Sunaev alpha parameter for disks", "/maps", unit=U.none, dependencies=[ "avg_map_rho_avg", "avg_map_T_mwavg", "avg_map_vr", "avg_map_vphi", ], ), "alpha_grav": Rule( self, self._alpha_grav, "Map of the graviational contrib to\ Shakura&Sunaev alpha parameter for disks", "/maps", unit=U.none, dependencies=["avg_map_coldens", "avg_map_T_mwavg"], ), "rho": Rule( self, self._rho, "Density slice", "/maps", unit=self.info["unit_density"], ), "speed_h": Rule( self, self._speed_h, "Horizontal speed slice wrt the line of sight", "/maps", unit=self.info["unit_velocity"], ), "speed_v": Rule( self, self._speed_v, "Vertical speed slice wrt the line of sight", "/maps", unit=self.info["unit_velocity"], ), "B_h": Rule( self, self._speed_h, "Horizontal slice of the magnetic field wrt the line of sight", "/maps", unit=self.info["unit_mag"], ), "B_v": Rule( self, self._speed_v, "Vertical slice of the magnetic field wrt the line of sight", "/maps", unit=self.info["unit_mag"], ), "T": Rule( self, self._temperature, "Temperature slice", "/maps", dependencies=["rho"], unit=self.info["unit_pressure"] / self.info["unit_density"], ), "levels": Rule( self, self._levels, "Max level within line of sight", "/maps" ), "jeans": Rule( self, self._jeans, "Jeans lenght slice", "/maps", dependencies=["rho", "T"], ), "jeans_ratio": Rule( self, self._jeans_ratio, "Jeans' lenght divided by the max resolution", "/maps", dependencies=["jeans", "levels"], ), "Q": Rule( self, self._toomreQ_disk, "Toomre Q parameter for a Keplerian disk", "/maps", dependencies=["coldens"], ), "sinks": Rule( self, self._sinks, group="/datasets", unit={ "Id": U.none, "M": U.Msun, "dmf": U.Msun, "x": "", "y": "", "z": "", "vx": "", "vy": "", "vz": "", "rot_period": "[y]", "lx": "|l|", "ly": "|l|", "lz": "|l|", "acc_rate": "[Msol/y]", "acc_lum": "[Lsol]", "age": U.year, "int_lum": "[Lsol]", "Teff": U.K, }, ), "pspec": Rule(self, self._pspec, "Power spectrum", "/hdf5"), "particles": Rule(self, self._particles, "Particles file", "/hdf5"), "filaments": Rule( self, self._filaments, "Filaments", "/datasets", dependencies={self.pp_params.filaments.datamap: "z"}, ), "filaments_forces": Rule( self, self._filaments_forces, "Filaments", "/datasets", dependencies={ "filaments": None, "slice_velphi": "z", "slice_gr": "z", "slice_P": "z", "slice_rho": "z", "slice_velr": "z", }, ), # Helpers "radial_bins": Rule( self, self._radial_bins, "Radial bins", "/radial", ), "radial_centers": Rule( self, self._radial_centers, "Centers of radial bins", "/radial", dependencies=["radial_bins"], ), "rr": Rule(self, self._rr, "Coordinate map", "/maps"), "bins_on_map": Rule( self, self._bins_on_map, "Convert map coordinates to bins", "/maps", dependencies=["radial_bins", "rr"], ), "B_int": Rule( self, self._B_int, "Magnetic intensity slice", "/maps", dependencies=["rho"], unit=self.info["unit_mag"], ), # PDF "rho_pdf": Rule( self, partial(self._vol_pdf, partial(simple_getter, "rho"), logbins=True), "Global rho-PDF", "/hist", unit=self.info["unit_density"], ), "rho_pdf_mw": Rule( self, partial( self._vol_pdf, partial(simple_getter, "rho"), weight_func=mass_func, logbins=True, ), "Global rho-PDF", "/hist", unit=self.info["unit_density"], ), "T_pdf": Rule( self, 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), "Global cos fluctuation-PDF", "/hist", dependencies=["mwa_speed"], unit=U.none, ), "Brho": Rule( self, self._Brho, "Average of B as a function of rho", "/datasets", unit={"rho": self.info["unit_density"], "B": self.info["unit_mag"]}, ), "Ek_Eb_rho": Rule( self, self._Ek_Eb_rho, "Average of Ek/Eb as a function of rho", "/datasets", dependencies=["mwa_speed"], unit={"rho": self.info["unit_density"], "Ek_Eb_rho": U.none}, ), # Profiles "axis": Rule( self, partial(self._get_axis), "Axis", "/profile", unit=self.info["unit_length"], ), "rho_prof": Rule( self, partial(self._plane_avg_uniform, partial(simple_getter, "rho")), "Rho profile", "/profile", unit=self.info["unit_density"], dependencies=["axis"], ), # globals "time_num": Rule( self, lambda: self.info["time"], "Time", "/globals", unit=self.info["unit_time"], ), "mass": Rule( self, partial(self._sum, mass_func), "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")), "Mass weighted speed average", "/globals", unit=self.info["unit_velocity"], ), "mwa_sigma": Rule( self, self._mwa_sigma, "Mass weighted speed average", "/globals", dependencies={"mwa_speed": None}, unit=self.info["unit_velocity"], ), "mwa_B_int": Rule( self, partial(self._vol_avg, getter_B_int), "Mass weighted Magnetic intensity average", "/globals", unit=self.info["unit_mag"], ), } # Average and other averageables = [ "coldens", "rho", "T", "Q", "vphi", "vr", "rho_avg", "P_avg", "T_avg", "P_mwavg", "T_mwavg", "alpha_disk", "alpha_grav", ] for name in averageables: self.rules["rad_avg_" + name] = Rule( self, partial(self._rad_avg, name), "Azimuthal average of {}".format(name), "/radial", dependencies=["radial_centers", "bins_on_map", name], ) self.rules["avg_map_" + name] = Rule( self, partial(self._rad_avg_map, name), "Interpolated map of azimuthal average of {}".format(name), "/maps", dependencies=["radial_centers", "bins_on_map", "rad_avg_" + name], ) self.rules["fluct_" + name] = Rule( self, partial(self._fluct_map, name), "Fluctuation wrt to average of {}".format(name), "/maps", dependencies=[name, "avg_map_" + name], ) self.rules["pdf_" + name] = Rule( self, partial(self._rad_fluct_pdf, name), "Probability density function of {} fluctuations".format(name), "/hist", dependencies=["rr", "fluct_" + name], ) self.rules["fit_pdf_" + name] = Rule( self, partial(self._fit_pdf, name), "Fit the PDF of {} fluctuations".format(name), "/hist", dependencies=["pdf_" + name], ) for name_bin in averageables: if name_bin is not name: self.rules["mbb_" + name + "_" + name_bin] = Rule( self, partial(self._mean_by_bins, name_bin, name), "Mean of {} by bins of {}".format(name, name_bin), "/hist", dependencies=[name, name_bin], ) self._gen_rule_transform("fluct_coldens", np.nanmax, "max", group="/globals") # Generic rules directly from Ramses fields for field in self.pp_params.pymses.variables: def generic_rule(name, getter, unit, oct_getter=None): if oct_getter is None: oct_getter = getter self.rules["slice_" + name] = Rule( self, partial(self._slice, getter, z=0.0, unit=unit), "{} slice".format(name), "/maps", unit=unit, ) self.rules[name + "_mwavg"] = Rule( self, partial(self._ax_avg, oct_getter, mass_weighted=True, unit=unit), "Ax mass-weighted averaged {}".format(name), "/maps", unit=unit, ) self.rules[name + "_avg"] = Rule( self, partial(self._ax_avg, oct_getter, mass_weighted=False, unit=unit), "Ax averaged {}".format(name), "/maps", unit=unit, ) # special for vectors if field in ["g", "vel"]: # Components for i, dir in enumerate(["x", "y", "z"]): generic_rule( field + dir, partial(vect_getter, field, i), self.unit_key[field], ) # Radial generic_rule( field + "r", partial(self.getter_vect_r, name_vect=field), self.unit_key[field], oct_getter=self.oct_getter_vect_r, ) # Othoradial generic_rule( field + "phi", partial(self.getter_vect_phi, name_vect=field), self.unit_key[field], oct_getter=self.oct_getter_vect_phi, ) # Norm generic_rule( field + "_norm", partial(norm_getter, field), self.unit_key[field] ) else: generic_rule(field, partial(simple_getter, field), self.unit_key[field]) super(PostProcessor, self).def_rules() def get_time(path, num): """ Return the time of the output (code units) Parameters ---------- num output number path_out path of the pipeline output Returns ------- time the time of the output (code units) """ try: f = open( path + "/output_" + str(num).zfill(5) + "/info_" + str(num).zfill(5) + ".txt" ) for line in f: ls = line.split() if len(ls) > 1 and ls[0] == "time": time = float(ls[2]) break f.close() return time except IOError as e: print(e) return np.nan