From beb6203d060bb6249face6840c606f69e5f85ab7 Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Thu, 16 Jan 2020 17:35:13 +0100 Subject: [PATCH] Possible to call rule directly --- .gitignore | 1 - plotter.py | 116 +++++++++++++--- postprocessor.py | 341 ++++++++++++++++++++++++++++++++++------------- pp_params.yml | 68 +++++----- pspec_new.py | 190 ++++++++++++++------------ run_selector.py | 250 ++++++++++++++++++---------------- 6 files changed, 626 insertions(+), 340 deletions(-) delete mode 100644 .gitignore diff --git a/.gitignore b/.gitignore deleted file mode 100644 index d5e992d..0000000 --- a/.gitignore +++ /dev/null @@ -1 +0,0 @@ -/.ipynb_checkpoints/ diff --git a/plotter.py b/plotter.py index d4d6dd6..4227b23 100644 --- a/plotter.py +++ b/plotter.py @@ -55,6 +55,7 @@ class Plotter(Aggregator, BaseProcessor): label_convert = { "turb_rms": "$f_{rms}$", "beta": "$\\beta$", + "beta_cool": "$\\beta_{c}$", "dens0": "$n_0$", "sfr_avg_window": "window", "comp_frac": "$\\zeta$", @@ -286,6 +287,7 @@ class Plotter(Aggregator, BaseProcessor): title=None, nml_key=None, put_time=True, + time_unit=cst.Myr, cmap="plasma", norm="log", autoscale=True, @@ -330,7 +332,9 @@ class Plotter(Aggregator, BaseProcessor): if put_time: time = self.save.root._v_attrs.time * self.comp.info["unit_time"] - time_str = "time = {:.6g} Myr".format(time.express(cst.Myr)) + time_str = "time = {:.6g} {}".format( + time.express(time_unit), time_unit.latex + ) if len(title) > 0: title = title + " | " + time_str else: @@ -423,19 +427,51 @@ class Plotter(Aggregator, BaseProcessor): if not label is None: P.ylabel(label) - def _plot_hist(self, name, ax_los="z", label=None, ylog=False): + def _plot_hist( + self, + name, + ax_los, + run, + label=None, + unit=None, + unit_coeff=1.0, + title=None, + nml_key=None, + put_time=True, + time_unit=cst.Myr, + ylog=False, + **kwargs + ): - pdf = self.save.get_node("/hist/" + name + "_" + ax_los) - values, centers = pdf.read() + node = self.save.get_node("/hist/" + name + "_" + ax_los) + + label, unit_old, unit = self._ax_label_unit(node, label, unit, unit_coeff) + values, centers = node.read() * unit_old.express(unit) width = centers[1] - centers[0] - P.bar(centers, values, width, log=ylog) + + P.bar(centers, values, width, log=ylog, **kwargs) P.grid() + + title = self._label_run(run, node, title, nml_key) + + if put_time: + time = self.save.root._v_attrs.time * self.comp.info["unit_time"] + time_str = "time = {:.6g} {}".format( + time.express(time_unit), time_unit.latex + ) + if len(title) > 0: + title = title + " | " + time_str + else: + title = time_str + + P.title(title) + if not label is None: P.xlabel(label) if "/hist/fit_" + name + "_" + ax_los in self.save: - slope = pdf.attrs.slope - origin = pdf.attrs.origin + slope = node.attrs.slope + origin = node.attrs.origin P.plot( centers, 10 ** (slope * centers + origin), @@ -481,6 +517,8 @@ class Plotter(Aggregator, BaseProcessor): if node_y._v_attrs.CLASS == "ARRAY": x = node_x.read() * xunit_old.express(xunit) y = node_y.read() * yunit_old.express(yunit) + mask = np.isfinite(y) + x, y = x[mask], y[mask] if smooth > 0: y = gaussian_filter1d(y, sigma=smooth) yerr = None @@ -488,6 +526,8 @@ class Plotter(Aggregator, BaseProcessor): elif "mean" in node_y: x = node_x.read() * xunit_old.express(xunit) y = node_y.mean.read() * yunit_old.express(yunit) + mask = np.isfinite(y) + x, y = x[mask], y[mask] if smooth > 0: y = gaussian_filter1d(y, sigma=smooth) yerr = node_y.std.read() * yunit_old.express(yunit) @@ -500,6 +540,8 @@ class Plotter(Aggregator, BaseProcessor): x_run, y_run = node_x[run], node_y[run] x = x_run.read() * xunit_old.express(xunit) y = y_run.read() * yunit_old.express(yunit) + mask = np.isfinite(y) + x, y = x[mask], y[mask] if smooth > 0: y = gaussian_filter1d(y, sigma=smooth) label_run = self._label_run(run, y_run, label, nml_key) @@ -507,16 +549,25 @@ class Plotter(Aggregator, BaseProcessor): P.legend() if linearfit: - _overlay_linearfit(x, y, yerr) + self._overlay_linearfit(x, y, yerr) - def _overlay_linearfit(x, y, yerr=None): + def _overlay_linearfit(self, x, y, yerr=None, fit_order=1): if yerr is None: - (a, b, rho, _, stderr) = linregress(x, y) + (a, b, rho, _map_rule, stderr) = linregress(x, y) + self._log( + "Linear fit y = {} x + {} with R^2 = {} and error is {}".format( + a, b, rho, stderr + ) + ) else: - c = polyfit(x, y, 1, w=[1.0 / ty for ty in yerr], full=False) + fit = polyfit(x, y, 1, w=[1.0 / ty for ty in yerr], full=True) + c = fit[0] + residual = fit[1][0][0] b, a = c[0], c[1] - - P.plot(x, a * y + b, "--", linewidth=1.5) + self._log( + "Linear fit y = {} x + {} with residual {}".format(a, b, residual) + ) + P.plot(x, a * x + b, "--", linewidth=1.5) def overlay_kennicutt(self, n0, step): P.grid(False) @@ -638,7 +689,7 @@ class Plotter(Aggregator, BaseProcessor): label="{}/avg({})".format(name, name), ), "Probability density function of {} fluctuations".format(name), - dependencies=["pdf_" + name], + dependencies=["fit_pdf_" + name], ) self.rules.update( @@ -647,12 +698,14 @@ class Plotter(Aggregator, BaseProcessor): self, partial( self._plot, - "/comp/beta", - "/comp/avg_pdf_slope_coldens", - linearfit=True, + "/comp/nml_cloud_params/beta_cool", + "/comp/avg_time_pdf_slope_coldens", ), kind="comp", - dependencies=["beta", "avg_pdf_slope_coldens"], + dependencies={ + "nml": "cloud_params/beta_cool", + "avg_time_pdf_slope_coldens": None, + }, ), "sink_mass": PlotRule( self, @@ -703,6 +756,17 @@ class Plotter(Aggregator, BaseProcessor): kind="series", dependencies=["rms_from_log"], ), + "turb_energy": PlotRule( + self, + partial( + self._plot, + "/series/rms_from_log/time", + "/series/rms_from_log/turb_energy", + xunit=cst.Myr, + ), + kind="series", + dependencies=["rms_from_log"], + ), "sigma": PlotRule( self, partial( @@ -713,15 +777,29 @@ class Plotter(Aggregator, BaseProcessor): xunit=cst.Myr, yunit=cst.km_s, ), - kind="comp", + kind="series", dependencies=["time_sigma"], ), + "max_fluct_coldens": PlotRule( + self, + partial( + self._plot, + "/series/time", + "/series/time_max_fluct_coldens_z", + ylabel="$\\max(\Sigma/\overline{\Sigma})$", + xunit=cst.Myr, + ), + kind="series", + dependencies={"time_max_fluct_coldens": "z"}, + ), "plot": PlotRule( self, lambda arg, **kwargs: self._plot(*arg, **kwargs), kind="comp" ), } ) + super(Plotter, self).def_rules() + class InteractiveGUI: """ diff --git a/postprocessor.py b/postprocessor.py index 6c159a5..20105f0 100644 --- a/postprocessor.py +++ b/postprocessor.py @@ -16,6 +16,9 @@ from pymses.sources.hop.file_formats import * from pymses.analysis import Camera, raytracing, slicing, splatting from pymses.filters import CellsToPoints from pymses.analysis import ScalarOperator, FractionOperator, MaxLevelOperator +import subprocess + +import module_extract as me from mypool import MyPool as Pool from functools import partial @@ -54,18 +57,17 @@ class Rule: return self.process_fn(**kwargs) def is_valid(self, arg): - save = self.postproc.save - valid = True - for dep in self.dependencies: - if dep in self.postproc.rules: - rule_dep = self.postproc.rules[dep] - if not arg is None: - valid = ( - valid and rule_dep.group + "/" + dep + "_" + str(arg) in save - ) - else: - valid = valid and rule_dep.group + "/" + dep in save - return valid and self.is_valid_add(arg) + # save = self.postproc.save + # valid = True + # for dep in self.dependencies: + # if dep in self.postproc.rules: + # rule_dep = self.postproc.rules[dep] + # if not arg is None: + # valid = valid and rule_dep.group + '/' + dep + '_' + str(arg) in save + # else: + # valid = valid and rule_dep.group + '/' + dep in save + # return valid and self.is_valid_add(arg) + return self.is_valid_add(arg) class BaseProcessor: @@ -142,17 +144,22 @@ class BaseProcessor: def _solve_and_process_rule(self, name, rule, arg, overwrite=False, **kwargs): updated = self._solve_dependencies(name, rule, arg, overwrite, **kwargs) overwrite_rule = overwrite or updated - self._process_rule(name, rule, arg, overwrite, **kwargs) + self._process_rule(name, rule, arg, overwrite_rule, **kwargs) def _solve_dependencies(self, name, rule, arg, overwrite=False, **kwargs): self.done_before_dep = len(self.just_done) - if type(rule.dependencies) == dict: - dep_arg = rule.dependencies[dep] - else: - dep_arg = arg + # Solve dependencies for dep in rule.dependencies: + try: + dep_arg = rule.dependencies[dep] + except: + dep_arg = arg + + if dep_arg == "__parent__": + dep_arg = arg + if self.solve_self_dep and dep in self.rules: rule_dep = self.rules[dep] self._solve_and_process_rule(dep, rule_dep, dep_arg, self.overwrite_dep) @@ -189,6 +196,53 @@ class BaseProcessor: else: self._log("{} is not valid in this context".format(name_full), "ERROR") + def def_rules(self): + for rule in self.rules: + setattr(self, rule, partial(self.process, rule)) + + +class HDF5Container(BaseProcessor): + + filename = "" + save = None + opened = False + + def open(self): + if not self.opened: + self.save = tables.open_file(self.filename, mode="a") + self.opened = True + + def close(self): + if self.opened: + self.save.close() + self.opened = False + + def _needs_computation(self, overwrite, name_full): + return overwrite or not (name_full in self.save) + + def _process_rule(self, name, rule, arg, overwrite, **kwargs): + self.open() + try: + super(HDF5Container, self)._process_rule( + name, rule, arg, overwrite, **kwargs + ) + finally: + self.close() + + def get_value(self, node_name): + self.open() + try: + node = self.save.get_node(node_name) + if node._v_attrs.CLASS == "GROUP": + value = {} + for child_name in node._v_children: + value[child_name] = self.get_value(node_name + "/" + child_name) + else: + value = node.read() + finally: + self.close() + return value + def _save_data(self, name_full, data, description, unit): """ Save data in the HDF5 structure, overwrite if necessary @@ -243,53 +297,13 @@ class BaseProcessor: ) self.save.get_node(name_full).attrs.unit = unit - @abstractmethod - def def_rules(self): - pass - - -class HDF5Container(BaseProcessor): - - filename = "" - save = None - opened = False - - def open(self): - if not self.opened: - self.save = tables.open_file(self.filename, mode="a") - self.opened = True - - def close(self): - if self.opened: - self.save.close() - self.opened = False - - def _needs_computation(self, overwrite, name_full): - return overwrite or not (name_full in self.save) - - def _process_rule(self, name, rule, arg, overwrite, **kwargs): + def set_value(self, node_name, data, description, unit): self.open() try: - super(HDF5Container, self)._process_rule( - name, rule, arg, overwrite, **kwargs - ) + self._save_data(node_name, data, description, unit) finally: self.close() - def get_value(self, node_name): - self.open() - try: - node = self.save.get_node(node_name) - if node._v_attrs.CLASS == "GROUP": - value = {} - for child_name in node._v_children: - value[child_name] = self.get_value(node_name + "/" + child_name) - else: - value = node.read() - finally: - self.close() - return value - def get_attribute(self, node_name, attr_name): self.open() try: @@ -301,9 +315,12 @@ class HDF5Container(BaseProcessor): def _map_rule(rule, arg, overwrite, path, path_out, pp_params, run_num): - pp = PostProcessor( - path + "/" + run_num[0], run_num[1], path_out + "/" + run_num[0], pp_params - ) + try: + pp = PostProcessor( + path + "/" + run_num[0], run_num[1], path_out + "/" + run_num[0], pp_params + ) + except pymses.RamsesIOError as e: + print(e) return pp.process(rule, arg, overwrite) @@ -395,9 +412,9 @@ class PostProcessor(HDF5Container): self._rt = raytracing.RayTracer(self._amr, self._ro.info, self._rho_op) # Set the extend of the image - self._radius = 0.5 / pp_params.out.zoom + self._radius = 0.5 / pp_params.pymses.zoom self._lbox = self.info["boxlen"] - center = pp_params.out.center + center = pp_params.pymses.center im_extent = [ (-self._radius + center[0]) * self._lbox, (self._radius + center[0]) * self._lbox, @@ -428,13 +445,13 @@ class PostProcessor(HDF5Container): ax_v = self._axes_v[ax_los] self._cam[ax_los] = Camera( - center=pp_params.out.center, + center=pp_params.pymses.center, line_of_sight_axis=ax_los, region_size=[2.0 * self._radius, 2.0 * self._radius], distance=self._radius, far_cut_depth=self._radius, up_vector=ax_v, - map_max_size=pp_params.out.map_size, + map_max_size=pp_params.pymses.map_size, ) self.close() @@ -673,7 +690,7 @@ class PostProcessor(HDF5Container): Computes the radius from the center """ im_extent = self.save.root.maps._v_attrs.im_extent - map_size = self.pp_params.out.map_size + map_size = self.pp_params.pymses.map_size pos_star = self.pp_params.disk.pos_star x = np.linspace(im_extent[0], im_extent[1], map_size) @@ -775,6 +792,31 @@ class PostProcessor(HDF5Container): pdf.attrs.var = np.var return True + def _clumps(self): + name = self.path_out + "/" + self.tag + "_" + str(self.num).zfill(5) + hop_save = name + "_hop" + "_prop_struct.save" + + me.make_clump_hop( + self.path, + self.num, + name + "_hop", + self.pp_params.hop.rho_thres, + self.pp_params.hop.lvl_thres, + [0.5, 0.5, 0.5], + 1, + path_out=path_out + "/", + path_hop="./", + force=True, + gcomp=False, + ) + hop_save = me.clump_properties( + name + "_hop", path, num, path_out=path_out + "/", gcomp=False + ) + f = open(path_out + "/" + hop_save) + hop_data = pickle.load(f) + f.close() + return hop_data + def _sinks(self): csv_name = ( self.path @@ -814,6 +856,56 @@ class PostProcessor(HDF5Container): return sinks_dict + def _transform(self, name, transform_fn, group="/maps", **kwargs): + src = self.save.get_node(group + "/" + name).read() + return transform_fn(src, **kwargs) + + def _gen_rule_transform( + self, + rule_src_name, + transform_fn, + transform_name, + subarray_name=None, + group=None, + ): + + rule_src = self.rules[rule_src_name] + + if subarray_name is None: + src_name = rule_src_name + group_src = rule_src.group + unit = rule_src.unit + description = rule_src.description + else: + src_name = subarray_name + group_src = rule_src.group + "/" + rule_src_name + unit = rule_src.unit[subarray_name] + description = rule_src.description[subarray_name] + + def fn(arg=None, **kwargs): + if arg is None: + return self._transform( + src_name, transform_fn, group=group_src, **kwargs + ) + else: + return self._transform( + src_name + "_" + str(arg), transform_fn, group=group_src, **kwargs + ) + + if group is None: + group = group_src + + name = transform_name + "_" + rule_src_name + + self.rules[name] = Rule( + self, + fn, + group=group, + unit=unit, + description=description, + dependencies=[rule_src_name], + ) + def def_rules(self): self.rules = { @@ -904,6 +996,7 @@ class PostProcessor(HDF5Container): "Teff": cst.K, }, ), + "clumps": Rule(self, self._clumps, group="/datasets"), # Helpers "radial_bins": Rule(self, self._radial_bins, "Radial bins", "/radial"), "rr": Rule(self, self._rr, "Coordinate map", "/maps"), @@ -980,6 +1073,10 @@ class PostProcessor(HDF5Container): dependencies=["pdf_" + name], ) + self._gen_rule_transform("fluct_coldens", np.max, "max", group="/globals") + + super(PostProcessor, self).def_rules() + class Comparator(Aggregator, HDF5Container): """ @@ -1033,6 +1130,7 @@ class Comparator(Aggregator, HDF5Container): ) self.namelist = selector.namelist + # Get info from one output. TODO Avoid using pymses for that self.info = self.pp_runs[self.runs[0]][self.nums[self.runs[0]][0]].info # log info @@ -1066,19 +1164,22 @@ class Comparator(Aggregator, HDF5Container): super(Comparator, self)._save_data(name_full, data, description, unit) self.save.get_node(name_full)._v_attrs.nums = self.nums - def _time_series(self, getter): + def _time_series(self, getter, arg=None): series = {} for run in self.runs: series[run] = np.zeros(len(self.nums[run])) for i, num in enumerate(self.nums[run]): - series[run][i] = getter(run, num) + series[run][i] = getter(run, num, arg=arg) return series - def _comp(self, getter): + def _comp(self, getter, use_num=True): prop = np.zeros(len(self.runs)) for i, run in enumerate(self.runs): - num = self.nums[run][0] - prop[i] = getter(run, num) + if use_num: + num = self.nums[run][0] + prop[i] = getter(run, num) + else: + prop[i] = getter(run) return prop def _time_avg(self, name, start=None, end=None, span=None, group="/series"): @@ -1108,21 +1209,22 @@ class Comparator(Aggregator, HDF5Container): std[i] = np.nanstd(serie[mask]) return {"runs": self.runs, "mean": mean, "std": std, "median": median} - def get_attr(self, attr_name, run, num, node_name="/"): + def get_attr(self, attr_name, run, num, node_name="/", arg=None): pp = self.pp_runs[run][num] + if not arg is None: + node_name = node_name + "_" + str(arg) return pp.get_attribute(node_name, attr_name) - def get_global(self, node_name, run, num, args=[None], unload_cells=False): + def get_global(self, node_name, run, num, arg=None, unload_cells=False): + if not arg is None: + node_name = node_name + "_" + str(arg) pp = self.pp_runs[run][num] if unload_cells: pp.unload_cells() value = pp.get_value(node_name) return value - def get_nml(self, nml_key, run, num=0): - """ - num parameter is ignored - """ + def get_nml(self, nml_key, run): res = self.namelist[run] for key in nml_key.split("/"): res = res[key] @@ -1135,7 +1237,7 @@ class Comparator(Aggregator, HDF5Container): return slope def _extract_sinks_from_log(self, series, log_filename, run): - cmd_grep = "grep 'Number of sink' {} -A 2".format(log_filename) + cmd_grep = "sed '/cpu.*/d' {} | grep 'Number of sink' -A 2".format(log_filename) content = os.popen(cmd_grep).readlines() for i in range(0, len(content), 4): series["nb_sink"][run].append(np.int(content[i].split("=")[1])) @@ -1154,11 +1256,18 @@ class Comparator(Aggregator, HDF5Container): return series def _extract_rms_from_log(self, series, log_filename, run): - cmd_grep = "grep 'turbulent rms' {} -B 1".format(log_filename) + cmd_grep = "grep 'turbulent rms' {} -C 1".format(log_filename) content = os.popen(cmd_grep).readlines() - for i in range(0, len(content), 3): + for i in range(0, len(content), 4): series["time"][run].append(np.float(content[i].split("=")[2].split()[0])) series["turb_rms"][run].append(np.float(content[i + 1].split(":")[1])) + try: + turb_energy = np.float(content[i + 2].split(":")[1]) + threshold = self.pp_params.rules.turb_energy_threshold + assert threshold < 0 or abs(turb_energy) < threshold + series["turb_energy"][run].append(turb_energy) + except (ValueError, IndexError, AssertionError): + series["turb_energy"][run].append(np.nan) return series def _from_log(self, keys, extractor): @@ -1223,6 +1332,33 @@ class Comparator(Aggregator, HDF5Container): ssfr[run][shift:] = sfr / surface return ssfr + def _gen_rule_time_global( + self, + glob_name, + name=None, + glob_group="/globals", + subarray_name=None, + unload_cells=True, + unit=cst.none, + description="", + ): + + if name is None: + name = "time_" + glob_name + + self.rules[name] = Rule( + self, + partial( + self._time_series, + partial( + self.get_global, glob_group + "/" + glob_name, unload_cells=True + ), + ), + group="/series", + unit=unit, + dependencies={"time": None, glob_name: "__parent__"}, + ) + def _gen_rule_avg(self, rule_src_name, subarray_name=None): rule_src = self.rules[rule_src_name] @@ -1307,11 +1443,21 @@ class Comparator(Aggregator, HDF5Container): "rms_from_log": Rule( self, partial( - self._from_log, ["time", "turb_rms"], self._extract_rms_from_log + self._from_log, + ["time", "turb_rms", "turb_energy"], + self._extract_rms_from_log, ), group="/series", - unit={"time": self.info["unit_time"], "turb_rms": cst.none}, - description={"time": "Time", "turb_rms": "Computed turbulent RMS"}, + unit={ + "time": self.info["unit_time"], + "turb_rms": cst.none, + "turb_energy": cst.none, + }, + description={ + "time": "Time", + "turb_rms": "Computed turbulent RMS", + "turb_energy": "Injected turbulent energy", + }, ), # Read from outputs "time": Rule( @@ -1323,31 +1469,46 @@ class Comparator(Aggregator, HDF5Container): unit=self.info["unit_time"], dependencies=["time_num"], ), - "time_sigma": Rule( + "time_pdf_slope_coldens": Rule( self, partial( self._time_series, - partial(self.get_global, "/globals/mwa_sigma", unload_cells=True), + partial( + self.get_attr, + "slope", + node_name="/hist/pdf_coldens_z", + ), ), group="/series", - unit=self.info["unit_velocity"], - dependencies=["time", "mwa_sigma"], + dependencies={"time": None, "fit_pdf_coldens": "z"}, ), # namelist "nml": Rule( self, - lambda nml_key: self._comp(partial(self.get_nml, nml_key)), + lambda nml_key: self._comp( + partial(self.get_nml, nml_key), use_num=False + ), group="/comp", ), } - for name in ["issfr", "time_sigma"]: + self._gen_rule_time_global( + "mwa_sigma", "time_sigma", unit=self.info["unit_velocity"] + ) + self._gen_rule_time_global("max_fluct_coldens") + + for name in ["issfr", "time_sigma", "time_pdf_slope_coldens"]: self._gen_rule_avg(name) self._gen_rule_avg("sinks_from_log", "mass_sink") self._gen_rule_avg("sinks_from_log", "nb_sink") self._gen_rule_avg("sfr_from_log", "sfr") + self._gen_rule_avg("rms_from_log", "turb_rms") + self._gen_rule_avg("rms_from_log", "turb_energy") + + super(Comparator, self).def_rules() + def get_time(path, num): """ diff --git a/pp_params.yml b/pp_params.yml index ef3c5ae..ab1b085 100644 --- a/pp_params.yml +++ b/pp_params.yml @@ -1,52 +1,58 @@ plot : # Plot parameters - out_ext : '.jpeg' # extension for plots - put_title : False # Add a title to plot - ntick : 6 # Number of ticks for maps - vel_red : 40 # Take point each vel_red for velocities + out_ext : '.jpeg' # extension for plots + put_title : False # Add a title to plot + # Maps + ntick : 6 # Number of ticks for maps + + # Overlays + vel_red : 40 # Take point each vel_red for velocities disk: # Disk speficic parameters - enable : False # Enable specific disk analysis - pos_star : [1., 1., 1.] # Position of the central star - binning : "log" # Kind of binning (lin : linear, log : logarithmic) - nb_bin : 100 # Number of bins for averaged quantities - bin_in : 1e-3 # Outer radius of the inner bin - bin_out : 0.25 # Inner radius of the outer bin - rmin_pdf : 0.075 # Inner radius for PDF computation - rmax_pdf : 0.3 # Outer radius for PDF computation + enable : False # Enable specific disk analysis + pos_star : [1., 1., 1.] # Position of the central star + binning : "log" # Kind of binning (lin : linear, log : logarithmic) + nb_bin : 100 # Number of bins for averaged quantities + bin_in : 1e-3 # Outer radius of the inner bin + bin_out : 0.25 # Inner radius of the outer bin + rmin_pdf : 0.075 # Inner radius for PDF computation + rmax_pdf : 0.3 # Outer radius for PDF computation pdf: # parameters for probability density functions - nb_bin : 50 # Number of bins for the PDF - xmin_fit : 0. # Lower boundary of the fit - xmax_fit : 1.25 # Upper boundary of the fit - + xmin_fit : 0. # Lower boundary of the fit + xmax_fit : 1.25 # Upper boundary of the fit pymses: # Parameters for Pymses reader - order : '<' # In which order the output are read - fft : False # Quick and dirty rendering using FFT - verbose : False # Let pymses write on standart output + # Source settings variables : ["rho","vel","Br","Bl","P", "g", "phi"] # Read these variables + order : '<' # Bit order + + # Processing options levelmax : 20 # Maximal AMR level visited when looking levels + fft : False # Quick and dirty rendering using FFT + verbose : False # Let pymses write on standart output -input: # Parameters on how to look for input files (: output from Ramses) - - log_prefix : "run.log" # Prefix of the log file - label_filename : "label.txt" # Name of the label file - nml_filename : "run.nml" # name of the namelist file - -out: # Parameters for post processing - + # Camera settings center : [0.5, 0.5, 0.5] # Center of the image zoom : 1. # Zoom of the image map_size : 1024 # Size of the computed maps in pixel - tag : "" # Tag for the image +input: # Parameters on how to look for input files (= output from Ramses) + log_prefix : "run.log" # Prefix of the log file + label_filename : "label.txt" # Name of the label file + nml_filename : "run.nml" # name of the namelist file + +out: # Parameters for post processing + tag : "" # Tag for the image interactive : False # Interactive mode (do not save the plots on the disk) -process: - verbose : True - num_process : 1 +process: # General setting of the post-processor module + verbose : True # Give more infos on what is going on + num_process : 1 # Number of forks + +rules: # Specific rules parameters + turb_energy_threshold : 1e13 # Remove invalid data (<0 = no threshold) diff --git a/pspec_new.py b/pspec_new.py index ed4cd11..4ec8536 100644 --- a/pspec_new.py +++ b/pspec_new.py @@ -10,16 +10,15 @@ import numpy as np from numpy.fft import fftn, ifft import pymses -from pymses.analysis import amr2cube +from pymses.analysis import ScalarOperator, Camera +from pymses.analysis import cube3d import pymses.utils.misc -# Don't use multiprocessing, it will crash with level 10 cubes -pymses.utils.misc.NUMBER_OF_PROCESSES_LIMIT = 1 - import tables as T +import bunch -from utils import args -from utils import tools +from i_utils import args +from i_utils import tools __generator__ = "pspec_new.py" __version__ = "0.2" @@ -420,85 +419,78 @@ def pspectrum(pcube, kcube, kbins, norm, nbinsf): return pspec, kbins, pspec2, fbins -if __name__ == "__main__": - # Command-line parser ---------------------------------------------------------- - parser = argparse.ArgumentParser( - formatter_class=argparse.RawDescriptionHelpFormatter, - description="Compute 2D and 3D power spectra", - epilog=textwrap.dedent( - """ +# Command-line parser ---------------------------------------------------------- +parser = argparse.ArgumentParser( + formatter_class=argparse.RawDescriptionHelpFormatter, + description="Compute 2D and 3D power spectra", + epilog=textwrap.dedent( + """ In the output file and node name formats, you can use the formatting fields: * %(iout)d - output number * %(varname)s - variable name * %(dim)d - dimension (3 for cube, 2 for slices) """ - ), - ) - parser.add_argument("repo", help="RAMSES output repository", type=str, default=".") - parser.add_argument( - "iouts", help="output numbers", type=args.selection, default=":" - ) - parser.add_argument( - "outfile", help="output file format (see below for fields)", type=str - ) - parser.add_argument( - "-n", - "--nodename", - help="node name format (see below for fields)", - type=str, - default="/out_%(iout)05d/d%(dim)d/%(varname)s", - ) - parser.add_argument( - "-O", - "--order", - help="byte order (= for native)", - type=str, - default="=", - choices=["<", ">", "="], - ) - parser.add_argument( - "-c", - "--center", - help="coordinates of the center", - type=args.center, - default=[0.5, 0.5, 0.5], - ) - parser.add_argument("-s", "--size", help="cube size", type=float, default=1.0) - parser.add_argument( - "-l", "--level", help="cube level (default: levelMIN)", type=int, default=0 - ) - parser.add_argument( - "-k", "--kbins", help="number of wave number bins", type=int, default=100 - ) - parser.add_argument( - "-f", - "--fbins", - help="number of Fourier bins for k-power 2D histogram (0 to disable)", - type=int, - default=0, - ) - parser.add_argument( - "-K", "--kbinsbig", help="number of big wave number bins", type=int, default=9 - ) - parser.add_argument( - "-d", - "--dkbig", - help="width of the big wave number bins", - type=float, - default=1.0, - ) - parser.add_argument( - "-S", - "--sliceaxis", - help="slicing axis", - type=str, - choices=["x", "y", "z"], - default="z", - ) + ), +) +parser.add_argument("repo", help="RAMSES output repository", type=str, default=".") +parser.add_argument("iouts", help="output numbers", type=args.selection, default=":") +parser.add_argument( + "outfile", help="output file format (see below for fields)", type=str +) +parser.add_argument( + "-n", + "--nodename", + help="node name format (see below for fields)", + type=str, + default="/out_%(iout)05d/d%(dim)d/%(varname)s", +) +parser.add_argument( + "-O", + "--order", + help="byte order (= for native)", + type=str, + default="=", + choices=["<", ">", "="], +) +parser.add_argument( + "-c", + "--center", + help="coordinates of the center", + type=args.center, + default=[0.5, 0.5, 0.5], +) +parser.add_argument("-s", "--size", help="cube size", type=float, default=1.0) +parser.add_argument( + "-l", "--level", help="cube level (default: levelMIN)", type=int, default=0 +) +parser.add_argument( + "-k", "--kbins", help="number of wave number bins", type=int, default=100 +) +parser.add_argument( + "-f", + "--fbins", + help="number of Fourier bins for k-power 2D histogram (0 to disable)", + type=int, + default=0, +) +parser.add_argument( + "-K", "--kbinsbig", help="number of big wave number bins", type=int, default=9 +) +parser.add_argument( + "-d", "--dkbig", help="width of the big wave number bins", type=float, default=1.0 +) +parser.add_argument( + "-S", + "--sliceaxis", + help="slicing axis", + type=str, + choices=["x", "y", "z"], + default="z", +) - arg = parser.parse_args() +def main(arg): add_pspec2 = False if arg.fbins > 0: add_pspec2 = True @@ -525,13 +517,11 @@ if __name__ == "__main__": if clvl == 0: clvl = ro.info["levelmin"] - read_lvl = max(clvl, ro.info["levelmin"]) + read_lvl = ro.info["levelmin"] # Extract cubes --------------------------------------------------------------- - xmin = [x - arg.size / 2.0 for x in arg.center] - xmax = [x + arg.size / 2.0 for x in arg.center] cube_vars = [ - "rho", + lambda dset: dset["rho"], lambda dset: dset["vel"][..., 0], lambda dset: dset["vel"][..., 1], lambda dset: dset["vel"][..., 2], @@ -539,11 +529,34 @@ if __name__ == "__main__": lambda dset: 0.5 * (dset["Bl"][..., 1] + dset["Br"][..., 1]), lambda dset: 0.5 * (dset["Bl"][..., 2] + dset["Br"][..., 2]), ] - cubes_arr = amr2cube(amr, cube_vars, xmin, xmax, read_lvl) + + cube_units = [ + ro.info["unit_density"], + ro.info["unit_velocity"], + ro.info["unit_velocity"], + ro.info["unit_velocity"], + ro.info["unit_mag"], + ro.info["unit_mag"], + ro.info["unit_mag"], + ] + + cam = Camera( + center=arg.center, + line_of_sight_axis="z", + region_size=[arg.size, arg.size], + distance=arg.size / 2.0, + far_cut_depth=arg.size / 2.0, + up_vector="y", + map_max_size=256, + ) + cubes = {} for i, v in enumerate(["rho", "vx", "vy", "vz", "Bx", "By", "Bz"]): - cubes[v] = cubes_arr[i, ...].copy() - del cubes_arr + operator = ScalarOperator(cube_vars[i], cube_units[i]) + extractor = cube3d.CubeExtractor(amr, operator) + cubes[v] = extractor.process( + cam, cube_size=arg.size, resolution=256 + ).data else: h5f = T.open_file("cube.hdf5", "r") read_lvl = np.asscalar(h5f.root.level.read()) @@ -836,3 +849,12 @@ if __name__ == "__main__": h5fd.create_array(metagrp, "version", __version__) h5fd.close() + + +if __name__ == "__main__": + arg = parser.parse_args() + main(arg) + + +def pspec(**kwargs): + main(bunch.bunchify(kwargs)) diff --git a/run_selector.py b/run_selector.py index 3e5f5f4..7dd27e9 100644 --- a/run_selector.py +++ b/run_selector.py @@ -7,138 +7,158 @@ from pp_params import * import f90nml - def __init__(self, path_in, - in_runs=None, in_nums="all", pp_params=default_params(), - name_run='*', namelist_cond = dict(), - sort_run_by=None, time_min=None, time_max=None): - self.path_in = path_in - self.pp_params = pp_params +class RunSelector: + def __init__( + self, + path_in, + in_runs=None, + in_nums="all", + pp_params=default_params(), + name_run="*", + namelist_cond={}, + sort_run_by=None, + time_min=None, + time_max=None, + ): + self.path_in = path_in + self.pp_params = pp_params - self.namelist = dict() - self.runs = self.get_runs(in_runs, name_run, namelist_cond, sort_run_by) + self.namelist = {} + self.runs = self.get_runs(in_runs, name_run, namelist_cond, sort_run_by) - self.info = dict() - for run in self.runs: - self.info[run] = dict() - self.nums = dict() + self.info = {} + for run in self.runs: + self.info[run] = {} + self.nums = {} - if not type(in_nums) == dict : - nums_temp = in_nums - in_nums = dict() - for run in self.runs: - in_nums[run] = nums_temp + if not type(in_nums) == dict: + nums_temp = in_nums + in_nums = {} + for run in self.runs: + in_nums[run] = nums_temp - for run in self.runs: - self.nums[run] = self.get_nums(run, in_nums[run], time_min, time_max) + for run in self.runs: + self.nums[run] = self.get_nums(run, in_nums[run], time_min, time_max) - def load_namelist(self, run): - path_run = self.path_in + "/" + run - path_nml = path_run + '/' + self.pp_params.input.nml_filename - return f90nml.read(path_nml) + def load_namelist(self, run): + path_run = self.path_in + "/" + run + path_nml = path_run + "/" + self.pp_params.input.nml_filename + return f90nml.read(path_nml) - def get_nml_value(self, nml_key, run): - res = self.namelist[run] - for key in nml_key.split('/'): - if key in res: - res = res[key] - elif key == nml_key.split('/')[-1]: - res = None - else: - raise KeyError(key) - return res + def get_nml_value(self, nml_key, run): + res = self.namelist[run] + for key in nml_key.split("/"): + if key in res: + res = res[key] + elif key == nml_key.split("/")[-1]: + res = None + else: + raise KeyError(key) + return res - def get_runs(self, in_runs=None, name_run='*', namelist_cond=dict(), sort_run_by=None): - def try_load_nml(run): - try : - self.namelist[run] = self.load_namelist(run) - success = True - except IOError: - success = False - return success + def get_runs(self, in_runs=None, name_run="*", namelist_cond={}, sort_run_by=None): + def try_load_nml(run): + try: + self.namelist[run] = self.load_namelist(run) + success = True + except IOError: + success = False + return success - runs = map(os.path.basename, filter(os.path.isdir, glob.glob(self.path_in + "/" + name_run))) - if not in_runs is None : - runs = filter(lambda n : n in runs, in_runs) - runs = filter(try_load_nml, runs) + runs = map( + os.path.basename, + filter(os.path.isdir, glob.glob(self.path_in + "/" + name_run)), + ) + if not in_runs is None: + runs = filter(lambda n: n in runs, in_runs) + runs = filter(try_load_nml, runs) - if type(namelist_cond) == tuple: - namelist_cond = [namelist_cond] + if type(namelist_cond) == tuple: + namelist_cond = [namelist_cond] - for (nml_key, operator, operand) in namelist_cond: - value = dict() - for run in runs: - value[run] = self.get_nml_value(nml_key, run) - if operator == '=': - runs = filter(lambda r: value[r] == operand, runs) - if operator == '!=': - runs = filter(lambda r: not value[r] == operand, runs) - elif operator == '>': - runs = filter(lambda r: value[r] > operand, runs) - elif operator == '<': - runs = filter(lambda r: value[r] < operand, runs) - elif operator == 'in': - runs = filter(lambda r: value[r] in operand, runs) + for (nml_key, operator, operand) in namelist_cond: + value = {} + for run in runs: + value[run] = self.get_nml_value(nml_key, run) + if operator == "=": + runs = filter(lambda r: value[r] == operand, runs) + if operator == "!=": + runs = filter(lambda r: not value[r] == operand, runs) + elif operator == ">": + runs = filter(lambda r: value[r] > operand, runs) + elif operator == "<": + runs = filter(lambda r: value[r] < operand, runs) + elif operator == "in": + runs = filter(lambda r: value[r] in operand, runs) + # Sort by the value in the namelist of sort_run_by + if not sort_run_by is None: + if type(sort_run_by) == str: + sort_run_by = [sort_run_by] + for nml_key in reversed(sort_run_by): + runs.sort(key=partial(self.get_nml_value, nml_key)) - # Sort by the value in the namelist of sort_run_by - if not sort_run_by is None: - if type(sort_run_by) == str: - sort_run_by = [sort_run_by] - for nml_key in reversed(sort_run_by): - runs.sort(key=partial(self.get_nml_value, nml_key)) + return runs - return runs + def load_info(self, run, num): + info_file = open( + self.path_in + + "/" + + run + + "/" + + "output_" + + str(num).zfill(5) + + "/" + + "info_" + + str(num).zfill(5) + + ".txt", + "r", + ) + info = {} + for line in info_file.readlines(): + parsed = yaml.safe_load(line.replace("=", ":")) + if type(parsed) == dict: + info.update(parsed) + info_file.close() + return info - def load_info(self, run, num): - info_file = open(self.path_in + "/" + run + "/" + - "output_" + str(num).zfill(5) + "/" + - "info_" + str(num).zfill(5) + ".txt", "r") - info = dict() - for line in info_file.readlines(): - parsed = yaml.safe_load(line.replace("=", ":")) - if type(parsed) == dict: - info.update(parsed) - info_file.close() - return info + def get_nums(self, run, in_nums=None, time_min=None, time_max=None): + def try_load_info(num): + try: + self.info[run][num] = self.load_info(run, num) + success = True + except IOError: + success = False + return success - def get_nums(self, run, in_nums=None, time_min=None, time_max=None): - def try_load_info(num): - try : - self.info[run][num] = self.load_info(run, num) - success = True - except IOError: - success = False - return success + names = glob.glob( + self.path_in + "/" + run + "/output_[0-9][0-9][0-9][0-9][0-9]" + ) + nums = map(lambda n: int(n.split("/")[-1].split("_")[1]), names) - names = glob.glob(self.path_in + "/" + run + "/output_[0-9][0-9][0-9][0-9][0-9]") - nums = map(lambda n : int(n.split('/')[-1].split('_')[1]), names) + if type(in_nums) == int: + in_nums = [in_nums] + if type(in_nums) == list: + nums = filter(lambda n: n in nums, in_nums) + nums = np.sort(nums) - if type(in_nums) == int: - in_nums = [in_nums] - if type(in_nums) == list: - nums = filter(lambda n : n in nums, in_nums) + if in_nums == "first": + i = 0 + while i < len(nums) - 1 and not try_load_info(nums[i]): + i = i + 1 + nums = [nums[i]] + elif in_nums == "last": + i = len(nums) - 1 + while i > 0 and not try_load_info(nums[i]): + i = i - 1 + nums = [nums[i]] + else: + nums = filter(try_load_info, nums) - nums = np.sort(nums) + if not time_min is None: + nums = filter(lambda n: self.info[run][n]["time"] >= time_min, nums) + if not time_max is None: + nums = filter(lambda n: self.info[run][n]["time"] <= time_max, nums) - if in_nums == "first": - i = 0 - while i < len(nums) - 1 and not try_load_info(nums[i]): - i = i + 1 - nums = [nums[i]] - elif in_nums == "last": - i = len(nums) - 1 - while i > 0 and not try_load_info(nums[i]): - i = i - 1 - nums = [nums[i]] - else: - nums = filter(try_load_info, nums) - - if not time_min is None: - nums = filter(lambda n: self.info[run][n]['time'] >= time_min, nums) - if not time_max is None: - nums = filter(lambda n: self.info[run][n]['time'] <= time_max, nums) - - - - return nums + return nums