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

Add possibility to launch an external rule

Browse Source 
+ some plotter improvements
+ some namelist improvements
+ add datacube extraction
This commit is contained in:
Noe Brucy
2022-08-23 10:10:42 +02:00
parent fcfd7a8dab
commit cdafe5cb61
5 changed files with 207 additions and 167 deletions
Download Patch File Download Diff File Expand all files Collapse all files
aggregator.py
+8 -4
View File
@@ -2,8 +2,10 @@
import numpy as np import numpy as np
from functools import partial from functools import partial
from baseprocessor import Rule
import snapshotprocessor import snapshotprocessor
try: try:
from mpi4py.futures import MPIPoolExecutor from mpi4py.futures import MPIPoolExecutor
@@ -25,9 +27,8 @@ def _map_aux(fun, path, path_out, params, run_num, **kwargs):
return fun(snap, **kwargs) return fun(snap, **kwargs)
def _map_rule(snap, rule, arg, overwrite, overwrite_dep): def _map_rule(snap, rule, **kwargs):
return snap.process(rule, arg, overwrite, overwrite_dep) return snap.process(rule, **kwargs)
class Aggregator: class Aggregator:
def get_snap_list(self, select=None): def get_snap_list(self, select=None):
@@ -41,8 +42,11 @@ class Aggregator:
def map(self, func, select=None, num_process=None, **kwargs): def map(self, func, select=None, num_process=None, **kwargs):
snaps = self.get_snap_list(select) if isinstance(func, Rule):
return self.map(_map_rule, select, num_process, rule=func, **kwargs)
snaps = self.get_snap_list(select)
if num_process is None: if num_process is None:
num_process = self.params.process.num_process num_process = self.params.process.num_process
baseprocessor.py
+56 -16
View File
@@ -14,19 +14,19 @@ from tables import HDF5ExtError
from params import default_params, load_params from params import default_params, load_params
from units import U from units import U
class Rule: class Rule:
def __init__( def __init__(
self, self,
postproc,
process, process,
description="", description="",
group="", group="",
dependencies=[], dependencies=[],
kind="snapshot", kind="snapshot",
unit=U.none, unit=U.none,
name="",
): ):
self.postproc = postproc self.name=name
self.process_fn = process self.process_fn = process
self.dependencies = dependencies self.dependencies = dependencies
self.group = group self.group = group
@@ -39,8 +39,6 @@ class Rule:
return self.process_fn(arg, **kwargs) return self.process_fn(arg, **kwargs)
else: else:
return self.process_fn(**kwargs) return self.process_fn(**kwargs)
class BaseProcessor: class BaseProcessor:
""" """
Base class for processors, should not be instanciated Base class for processors, should not be instanciated
@@ -82,20 +80,39 @@ class BaseProcessor:
arg=None, arg=None,
overwrite=False, overwrite=False,
overwrite_dep=False, overwrite_dep=False,
skip_dep=False,
select=None, select=None,
**kwargs, **kwargs,
): ):
"""
Process the rule `to_process`
"""
self.overwrite_dep = overwrite_dep self.overwrite_dep = overwrite_dep
self.just_done = [] self.just_done = []
""" Process the rule 'to_process'
Parameters
----------
to_process : str of Rule
name of the rule to process or Rule object with nonempty rule.name
arg : optional
argument to give to the rule
overwrite : bool, optional
Force redo if already done
overwrite_dep : bool, optional
Force redoing of the dependencies even if already done
skip_dep : bool, optional
Skip the dependency checks (assume they are already done)
select : dict, optional
Select object (see RunSelector) to only select some run/snapshot
"""
if to_process in self.rules: if to_process in self.rules:
rule = self.rules[to_process] rule = self.rules[to_process]
return self._solve_and_process_rule( return self._solve_and_process_rule(
to_process, rule, arg, overwrite, select, **kwargs to_process, rule, arg, overwrite, skip_dep, select, **kwargs
)
elif isinstance(to_process, Rule):
rule = to_process
return self._solve_and_process_rule(
rule.name, rule, arg, overwrite, skip_dep, select, **kwargs
) )
else: else:
self._log( self._log(
@@ -106,9 +123,30 @@ class BaseProcessor:
) )
def _solve_and_process_rule( def _solve_and_process_rule(
self, name, rule, arg, overwrite=False, select=None, **kwargs self, name, rule, arg, overwrite=False, skip_dep=False, select=None, **kwargs
): ):
updated = self._solve_dependencies(name, rule, arg, overwrite, select) """Resolve dependencies and proceed in the processing of a rule
Parameters
----------
name : str
name of the rule
rule : Rule
rule object
overwrite : bool, optional
Force redo if already done
skip_dep : bool, optional
Skip the dependency checks (assume they are already done)
select : dict, optional
Select object (see RunSelector) to only select some run/snapshot
Returns
-------
The outbut of self._process_rule
"""
updated = False
if not skip_dep:
updated = self._solve_dependencies(name, rule, arg, overwrite, select)
overwrite_rule = overwrite or updated overwrite_rule = overwrite or updated
return self._process_rule(name, rule, arg, overwrite_rule, select, **kwargs) return self._process_rule(name, rule, arg, overwrite_rule, select, **kwargs)
@@ -147,7 +185,7 @@ class BaseProcessor:
return overwrite return overwrite
def _process_rule(self, name, rule, arg, overwrite=False, select=None, **kwargs): def _process_rule(self, name, rule, arg, overwrite=False, select=None, **kwargs):
if arg is not None: if arg is not None and not isinstance(arg, BaseProcessor):
name_full = rule.group + "/" + name + "_" + str(arg) name_full = rule.group + "/" + name + "_" + str(arg)
else: else:
name_full = rule.group + "/" + name name_full = rule.group + "/" + name
@@ -418,7 +456,6 @@ class HDF5Container(BaseProcessor):
name = transform_name + "_" + rule_src_name name = transform_name + "_" + rule_src_name
self.rules[name] = Rule( self.rules[name] = Rule(
self,
fn, fn,
group=group, group=group,
unit=unit, unit=unit,
@@ -468,4 +505,7 @@ def oct_vect_getter(name, i, dset):
def norm_getter(name, dset): def norm_getter(name, dset):
return np.sqrt(np.sum(dset[name] ** 2, axis=1)) return np.sqrt(np.sum(dset[name] ** 2, axis=1))
def oct_norm_getter(name, dset):
return np.sqrt(np.sum(dset[name] ** 2, axis=2))
plotter.py
+42 -64
View File
@@ -155,8 +155,7 @@ def line_integral_convolution(ax, map_h, map_v, extent, **kwargs):
class PlotRule(Rule): class PlotRule(Rule):
""" """
The rule class, speficic to plot. The rule class, specific to plot.
Add an extra method, plot, that take the reference to an open hdf5 file (from pytables)
""" """
def datafile(self, name, arg): def datafile(self, name, arg):
@@ -197,10 +196,16 @@ class Plotter(Aggregator, BaseProcessor):
"comp_frac": "$\chi$", "comp_frac": "$\chi$",
} }
# Conversion table from namelist values (from amses config file) into LaTex strings # Conversion table from namelist values (from ramses config file) into LaTex strings
value_convert = { value_str = {
"sfr_avg_window": lambda x: "${:g}$ Myr".format(80 * x), "sfr_avg_window": lambda x: "${:g}$ Myr".format(80 * x),
"Bx": lambda x: "${:g}$ $\\mu G$".format(7.6189439 * x), "Bx": lambda x: "${:.1f}$ $\\mu G$".format(7.6189439 * x),
}
# Conversion table from namelist values (from ramses config file) into suitanle units
value_convert = {
"sfr_avg_window": lambda x: 80 * x, # Myr
"Bx": lambda x: x * 7.6189439,
} }
def __init__( def __init__(
@@ -560,11 +565,16 @@ class Plotter(Aggregator, BaseProcessor):
prop_label = self.label_convert[prop_name] prop_label = self.label_convert[prop_name]
else: else:
prop_label = prop_name prop_label = prop_name
prop_value = self.study.get_nml(nml_key, run) try:
if prop_name in self.value_convert: prop_value = self.study.get_nml(nml_key, run)
except KeyError:
return ""
if prop_name in self.value_str:
prop_value_str = self.value_str[prop_name](prop_value)
elif prop_name in self.value_convert:
prop_value_str = self.value_convert[prop_name](prop_value) prop_value_str = self.value_convert[prop_name](prop_value)
elif type(prop_value) in [int, float]: elif type(prop_value) in [int, float]:
prop_value_str = convert_exp(prop_value, digits=5) prop_value_str = convert_exp(prop_value, digits=4)
else: else:
prop_value_str = str(prop_value) prop_value_str = str(prop_value)
return r"{} = {}".format(prop_label, prop_value_str) return r"{} = {}".format(prop_label, prop_value_str)
@@ -584,7 +594,9 @@ class Plotter(Aggregator, BaseProcessor):
elif nml_key is not None: elif nml_key is not None:
if not type(nml_key) == list: if not type(nml_key) == list:
nml_key = [nml_key] nml_key = [nml_key]
label_run = ", ".join(map(get_label_nml, nml_key)) lbl_list = map(get_label_nml, nml_key) # get namelist value
lbl_list = filter(lambda x: len(x) > 0, lbl_list) # Remove void labels
label_run = ", ".join(lbl_list)
if label is not None: if label is not None:
label_run = label + " (" + label_run + ")" label_run = label + " (" + label_run + ")"
@@ -667,6 +679,7 @@ class Plotter(Aggregator, BaseProcessor):
colorbar=True, colorbar=True,
embeded=False, embeded=False,
text_embeded=None, text_embeded=None,
text_kwargs={},
colorbar_embeded=None, colorbar_embeded=None,
axes_indicator=None, axes_indicator=None,
overtext_color="w", overtext_color="w",
@@ -814,10 +827,11 @@ class Plotter(Aggregator, BaseProcessor):
if text_embeded: if text_embeded:
ax.text( ax.text(
x=0.05, x=0.05,
y=0.95, y=0.91,
s=title, s=title,
color=overtext_color, color=overtext_color,
transform=ax.transAxes, transform=ax.transAxes,
**text_kwargs
) )
else: else:
plt.title(title) plt.title(title)
@@ -1021,7 +1035,7 @@ class Plotter(Aggregator, BaseProcessor):
c = c(data)[mask] c = c(data)[mask]
# Scatter plot # Scatter plot
plt.scatter(part_h, part_v, s=s, c=c, **kwargs) scatter = plt.scatter(part_h, part_v, s=s, c=c, **kwargs)
def _overlay_vector( def _overlay_vector(
self, self,
@@ -1143,11 +1157,13 @@ class Plotter(Aggregator, BaseProcessor):
).express(unit_time) ).express(unit_time)
color = colors(time) color = colors(time)
else: else:
nml = self.study.get_nml(nml_color, run) nml_value = self.study.get_nml(nml_color, run)
if os.path.basename(nml_color) in self.value_convert:
nml_value = self.value_convert[ os.path.basename(nml_color)](nml_value)
try: try:
color = colors[nml] color = colors[nml_value]
except TypeError: except TypeError:
color = colors(nml) color = colors(nml_value)
# Actual plot # Actual plot
if kind == "bar": if kind == "bar":
@@ -1209,6 +1225,7 @@ class Plotter(Aggregator, BaseProcessor):
label=None, label=None,
xunit=None, xunit=None,
yunit=None, yunit=None,
put_units=True,
xunit_coeff=1.0, xunit_coeff=1.0,
yunit_coeff=1.0, yunit_coeff=1.0,
xtransform=None, xtransform=None,
@@ -1255,10 +1272,10 @@ class Plotter(Aggregator, BaseProcessor):
# Find proper labels # Find proper labels
xlabel, xunit_old, xunit = self._ax_label_unit( xlabel, xunit_old, xunit = self._ax_label_unit(
name_x, xlabel, xunit, xunit_coeff name_x, xlabel, xunit, xunit_coeff, put_units=put_units,
) )
ylabel, yunit_old, yunit = self._ax_label_unit( ylabel, yunit_old, yunit = self._ax_label_unit(
name_y, ylabel, yunit, yunit_coeff name_y, ylabel, yunit, yunit_coeff, put_units=put_units,
) )
# If relevent, get time # If relevent, get time
@@ -1350,11 +1367,13 @@ class Plotter(Aggregator, BaseProcessor):
).express(unit_time) ).express(unit_time)
color = colors(time) color = colors(time)
else: else:
nml = self.study.get_nml(nml_color, run) nml_value = self.study.get_nml(nml_color, run)
if os.path.basename(nml_color) in self.value_convert:
nml_value = self.value_convert[os.path.basename(nml_color)](nml_value)
try: try:
color = colors[nml] color = colors[nml_value]
except TypeError: except TypeError:
color = colors(nml) color = colors(nml_value)
if yerr is None or yerr_kind is None: if yerr is None or yerr_kind is None:
(base_line,) = plt.plot(x, y, label=label, color=color, **kwargs) (base_line,) = plt.plot(x, y, label=label, color=color, **kwargs)
else: else:
@@ -1469,7 +1488,6 @@ class Plotter(Aggregator, BaseProcessor):
else: else:
name_rule = name_y + "_" + name_x name_rule = name_y + "_" + name_x
self.rules[name_rule] = PlotRule( self.rules[name_rule] = PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/series/" + logrule + "/" + name_x, "/series/" + logrule + "/" + name_x,
@@ -1485,20 +1503,15 @@ class Plotter(Aggregator, BaseProcessor):
This is where rules are defined This is where rules are defined
""" """
self.rules = { self.rules = {
"plot": PlotRule( "plot": PlotRule(lambda arg, **kwargs: self._plot(*arg, **kwargs), kind="comp"
self, lambda arg, **kwargs: self._plot(*arg, **kwargs), kind="comp"
), ),
"plot_run": PlotRule( "plot_run": PlotRule(lambda arg, **kwargs: self._plot(*arg, **kwargs), kind="run"
self, lambda arg, **kwargs: self._plot(*arg, **kwargs), kind="run"
), ),
"plot_snapshot": PlotRule( "plot_snapshot": PlotRule(lambda arg, **kwargs: self._plot(*arg, **kwargs)
self, lambda arg, **kwargs: self._plot(*arg, **kwargs)
), ),
"plot_map": PlotRule( "plot_map": PlotRule(lambda mapname, **kwargs: self._plot_map(mapname, **kwargs)
self, lambda mapname, **kwargs: self._plot_map(mapname, **kwargs)
), ),
"coldens": PlotRule( "coldens": PlotRule(
self,
partial( partial(
self._plot_map, self._plot_map,
"coldens", "coldens",
@@ -1509,7 +1522,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["coldens"], dependencies=["coldens"],
), ),
"slice_T": PlotRule( "slice_T": PlotRule(
self,
partial( partial(
self._plot_map, self._plot_map,
"T", "T",
@@ -1519,19 +1531,16 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["T"], dependencies=["T"],
), ),
"alpha_disk": PlotRule( "alpha_disk": PlotRule(
self,
partial(self._plot_map, "alpha_disk", label=r"$\alpha$"), partial(self._plot_map, "alpha_disk", label=r"$\alpha$"),
"Map of the Shakura&Sunaev alpha parameter for disks", "Map of the Shakura&Sunaev alpha parameter for disks",
dependencies=["alpha_disk"], dependencies=["alpha_disk"],
), ),
"alpha_grav": PlotRule( "alpha_grav": PlotRule(
self,
partial(self._plot_map, "alpha_grav", label=r"$\alpha_g$"), partial(self._plot_map, "alpha_grav", label=r"$\alpha_g$"),
"Map of the grav Shakura&Sunaev alpha parameter for disks", "Map of the grav Shakura&Sunaev alpha parameter for disks",
dependencies=["alpha_grav"], dependencies=["alpha_grav"],
), ),
"coldens_l": PlotRule( "coldens_l": PlotRule(
self,
partial( partial(
self._plot_map, self._plot_map,
"coldens", "coldens",
@@ -1543,7 +1552,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["coldens", "levels"], dependencies=["coldens", "levels"],
), ),
"slice_rho_v": PlotRule( "slice_rho_v": PlotRule(
self,
partial( partial(
self._plot_map, self._plot_map,
"slice_rho", "slice_rho",
@@ -1555,7 +1563,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["slice_rho"], dependencies=["slice_rho"],
), ),
"jeans_ratio": PlotRule( "jeans_ratio": PlotRule(
self,
partial( partial(
self._plot_map, self._plot_map,
"jeans_ratio", "jeans_ratio",
@@ -1568,7 +1575,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["jeans_ratio", "levels"], dependencies=["jeans_ratio", "levels"],
), ),
"Q": PlotRule( "Q": PlotRule(
self,
partial( partial(
self._plot_map, self._plot_map,
"Q", "Q",
@@ -1581,25 +1587,21 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["Q"], dependencies=["Q"],
), ),
"rho_pdf": PlotRule( "rho_pdf": PlotRule(
self,
partial(self._plot_hist, "rho_pdf"), partial(self._plot_hist, "rho_pdf"),
"$\rho$-PDF", "$\rho$-PDF",
dependencies=["rho_pdf"], dependencies=["rho_pdf"],
), ),
"rho_pdf_mw": PlotRule( "rho_pdf_mw": PlotRule(
self,
partial(self._plot_hist, "rho_pdf_mw"), partial(self._plot_hist, "rho_pdf_mw"),
"Mass weighted $\rho$-PDF", "Mass weighted $\rho$-PDF",
dependencies=["rho_pdf_mw"], dependencies=["rho_pdf_mw"],
), ),
"cos_pdf": PlotRule( "cos_pdf": PlotRule(
self,
partial(self._plot_hist, "cos_pdf"), partial(self._plot_hist, "cos_pdf"),
"cos-PDF", "cos-PDF",
dependencies=["cos_pdf", "mwa_speed"], dependencies=["cos_pdf", "mwa_speed"],
), ),
"avg_coldens_pdf": PlotRule( "avg_coldens_pdf": PlotRule(
self,
partial( partial(
self._plot_hist, self._plot_hist,
"avg_time_coldens_pdf_z", "avg_time_coldens_pdf_z",
@@ -1612,19 +1614,16 @@ class Plotter(Aggregator, BaseProcessor):
dependencies={"avg_time_coldens_pdf": "z"}, dependencies={"avg_time_coldens_pdf": "z"},
), ),
"T_pdf": PlotRule( "T_pdf": PlotRule(
self,
partial(self._plot_hist, "T_pdf"), partial(self._plot_hist, "T_pdf"),
"T-PDF on a 2D slice", "T-PDF on a 2D slice",
dependencies=["T_pdf"], dependencies=["T_pdf"],
), ),
"P_pdf": PlotRule( "P_pdf": PlotRule(
self,
partial(self._plot_hist, "P_pdf"), partial(self._plot_hist, "P_pdf"),
"P-PDF on a 2D slice ", "P-PDF on a 2D slice ",
dependencies=["P_pdf"], dependencies=["P_pdf"],
), ),
"Brho": PlotRule( "Brho": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/datasets/Brho/rho", "/datasets/Brho/rho",
@@ -1636,7 +1635,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["Brho"], dependencies=["Brho"],
), ),
"Ek_Eb_rho": PlotRule( "Ek_Eb_rho": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/datasets/Ek_Eb_rho/rho", "/datasets/Ek_Eb_rho/rho",
@@ -1648,14 +1646,12 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["Ek_Eb_rho", "mwa_speed"], dependencies=["Ek_Eb_rho", "mwa_speed"],
), ),
"rho_prof": PlotRule( "rho_prof": PlotRule(
self,
partial(self._plot, "/profile/axis", "/profile/rho_prof"), partial(self._plot, "/profile/axis", "/profile/rho_prof"),
"Density profile", "Density profile",
dependencies=["axis", "rho_prof"], dependencies=["axis", "rho_prof"],
), ),
"pspec": PlotRule(self, self._pspec, dependencies={"pspec": None}), "pspec": PlotRule(self, self._pspec, dependencies={"pspec": None}),
"sbeta": PlotRule( "sbeta": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/comp/nml_cloud_params/beta_cool", "/comp/nml_cloud_params/beta_cool",
@@ -1669,7 +1665,6 @@ class Plotter(Aggregator, BaseProcessor):
}, },
), ),
"sbeta_onavg": PlotRule( "sbeta_onavg": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/comp/sbeta_onavg/beta", "/comp/sbeta_onavg/beta",
@@ -1681,7 +1676,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["sbeta_onavg"], dependencies=["sbeta_onavg"],
), ),
"sink_mass": PlotRule( "sink_mass": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/series/sinks_from_log/time", "/series/sinks_from_log/time",
@@ -1694,7 +1688,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["sinks_from_log"], dependencies=["sinks_from_log"],
), ),
"ssm": PlotRule( "ssm": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/series/sinks_from_log/time", "/series/sinks_from_log/time",
@@ -1707,7 +1700,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["ssm"], dependencies=["ssm"],
), ),
"assfr": PlotRule( "assfr": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/series/sfr_from_log/time", "/series/sfr_from_log/time",
@@ -1720,7 +1712,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["sfr_from_log"], dependencies=["sfr_from_log"],
), ),
"issfr": PlotRule( "issfr": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/series/sinks_from_log/time", "/series/sinks_from_log/time",
@@ -1733,7 +1724,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["issfr"], dependencies=["issfr"],
), ),
"turb_rms": PlotRule( "turb_rms": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/series/rms_from_log/time", "/series/rms_from_log/time",
@@ -1745,7 +1735,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["rms_from_log"], dependencies=["rms_from_log"],
), ),
"turb_energy": PlotRule( "turb_energy": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/series/rms_from_log/time", "/series/rms_from_log/time",
@@ -1757,7 +1746,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["rms_from_log"], dependencies=["rms_from_log"],
), ),
"turb_power": PlotRule( "turb_power": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/series/rms_from_log/time", "/series/rms_from_log/time",
@@ -1769,7 +1757,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["turb_power"], dependencies=["turb_power"],
), ),
"sigma": PlotRule( "sigma": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/series/time", "/series/time",
@@ -1783,7 +1770,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["time_sigma"], dependencies=["time_sigma"],
), ),
"mwa_B_int": PlotRule( "mwa_B_int": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/series/time", "/series/time",
@@ -1796,7 +1782,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["time_mwa_B_int"], dependencies=["time_mwa_B_int"],
), ),
"mass": PlotRule( "mass": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/series/time", "/series/time",
@@ -1809,7 +1794,6 @@ class Plotter(Aggregator, BaseProcessor):
dependencies=["time_mass"], dependencies=["time_mass"],
), ),
"max_fluct_coldens": PlotRule( "max_fluct_coldens": PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/series/time", "/series/time",
@@ -1879,7 +1863,6 @@ class Plotter(Aggregator, BaseProcessor):
for name in averageables: for name in averageables:
self.rules["rad_" + name] = PlotRule( self.rules["rad_" + name] = PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/radial/radial_centers", "/radial/radial_centers",
@@ -1890,7 +1873,6 @@ class Plotter(Aggregator, BaseProcessor):
) )
self.rules["dispersion_rad_" + name] = PlotRule( self.rules["dispersion_rad_" + name] = PlotRule(
self,
partial( partial(
self._plot, self._plot,
"/radial/radial_centers", "/radial/radial_centers",
@@ -1901,28 +1883,24 @@ class Plotter(Aggregator, BaseProcessor):
) )
self.rules["avg_map_" + name] = PlotRule( self.rules["avg_map_" + name] = PlotRule(
self,
partial(self._plot_map, "avg_map_" + name), partial(self._plot_map, "avg_map_" + name),
"Map of the radial average of {}".format(name), "Map of the radial average of {}".format(name),
dependencies=["avg_map_" + name], dependencies=["avg_map_" + name],
) )
self.rules["mwavg_map_" + name] = PlotRule( self.rules["mwavg_map_" + name] = PlotRule(
self,
partial(self._plot_map, "mwavg_map_" + name), partial(self._plot_map, "mwavg_map_" + name),
"Map of the mass weighted radial average of {}".format(name), "Map of the mass weighted radial average of {}".format(name),
dependencies=["avg_map_" + name], dependencies=["avg_map_" + name],
) )
self.rules["fluct_" + name] = PlotRule( self.rules["fluct_" + name] = PlotRule(
self,
partial(self._plot_map, "fluct_" + name, cmap="RdBu_r"), partial(self._plot_map, "fluct_" + name, cmap="RdBu_r"),
"Fluctuation of {}".format(name), "Fluctuation of {}".format(name),
dependencies=["fluct_" + name], dependencies=["fluct_" + name],
) )
self.rules["pdf_" + name] = PlotRule( self.rules["pdf_" + name] = PlotRule(
self,
partial(self._plot_hist, "pdf_" + name, ylog=True), partial(self._plot_hist, "pdf_" + name, ylog=True),
"Probability density function of {} fluctuations".format(name), "Probability density function of {} fluctuations".format(name),
dependencies=["pdf_" + name], dependencies=["pdf_" + name],
snapshotprocessor.py
+69 -58
View File
@@ -27,6 +27,7 @@ from pymses.analysis import (
raytracing, raytracing,
slicing, slicing,
splatting, splatting,
cube3d,
) )
from pymses.filters import CellsToPoints, RegionFilter from pymses.filters import CellsToPoints, RegionFilter
from pymses.sources.hop.hop import HOP from pymses.sources.hop.hop import HOP
@@ -41,6 +42,7 @@ from baseprocessor import (
HDF5Container, HDF5Container,
Rule, Rule,
norm_getter, norm_getter,
oct_norm_getter,
simple_getter, simple_getter,
vect_getter, vect_getter,
oct_vect_getter, oct_vect_getter,
@@ -379,6 +381,14 @@ class SnapshotProcessor(HDF5Container):
else: else:
self.fil = None self.fil = None
# Convert time unit
if isinstance(unit_time, str):
factor = self.get_nml(unit_time)
unit_time = U.Unit(
name=os.path.basename(unit_time),
base_unit=self.info["unit_time"],
coeff=factor)
time_in_right_unit = self.time * self.info["unit_time"].express(unit_time) time_in_right_unit = self.time * self.info["unit_time"].express(unit_time)
self.snapshot = Snapshot( self.snapshot = Snapshot(
name=str(self.num), name=str(self.num),
@@ -636,6 +646,46 @@ class SnapshotProcessor(HDF5Container):
""" Azimuthal velocity """ """ Azimuthal velocity """
return self.oct_getter_vect_phi(dset, "vel") return self.oct_getter_vect_phi(dset, "vel")
def datacube(self, getter, level=None, unit=U.none):
"""
Return a datacube of the source box.
Parameters
----------
getter : callable
A callable that extract the wanted data from a pymses dataset
unit : U.Unit
Unit of the resulting dataset
Returns
-------
A numpy array containing the slice
"""
unit = self._get_units(unit)
operator = ScalarOperator(getter, unit)
extractor = cube3d.CubeExtractor(self._amr, operator)
if level is None:
level = self.get_nml("amr_params/levelmin")
size = 1.0
cam = Camera(
center=self.params.pymses.center,
line_of_sight_axis="z",
region_size=[size, size],
distance=size / 2.0,
far_cut_depth=size / 2.0,
up_vector="y",
map_max_size=2 ** level,
)
cube = extractor.process(
cam, cube_size=1.0, resolution=2 ** level
).data
return cube
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.
@@ -664,7 +714,8 @@ 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, oct_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
@@ -673,12 +724,12 @@ class SnapshotProcessor(HDF5Container):
""" """
unit = self._get_units(unit) unit = self._get_units(unit)
if surf_qty: if surf_qty:
op = ScalarOperator(getter, unit) op = ScalarOperator(oct_getter, unit)
else: else:
if mass_weighted: if mass_weighted:
def num(dset): def num(dset):
value = getter(dset) value = oct_getter(dset)
rho = getter_rho(dset) rho = getter_rho(dset)
return rho * value return rho * value
@@ -688,7 +739,7 @@ class SnapshotProcessor(HDF5Container):
else: # Volume weighted else: # Volume weighted
def num(dset): def num(dset):
value = getter(dset) value = oct_getter(dset)
return value return value
def denum(dset): def denum(dset):
@@ -779,9 +830,11 @@ 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, old_unit=None, unit=None, logbins=False, weight_func=vol_func):
self.load_cells() self.load_cells()
data = getter(self.cells) data = getter(self.cells)
if old_unit is not None and unit is not None:
data *= old_unit.express(unit)
if logbins: if logbins:
data = np.log10(data) data = np.log10(data)
weights = weight_func(self.cells) weights = weight_func(self.cells)
@@ -987,7 +1040,7 @@ class SnapshotProcessor(HDF5Container):
dmap_omega = rt_omega.process(self._cam[ax_los]).map.T dmap_omega = rt_omega.process(self._cam[ax_los]).map.T
return dmap_omega return dmap_omega
def _toomreQ_disk(self, ax_los, omega_approx=True, G1_units=True, coarsen_factor=1): def _toomreQ_disk(self, ax_los, omega_approx=False, G1_units=False, coarsen_factor=1):
""" """
Compute the Toomre Q parameter Compute the Toomre Q parameter
""" """
@@ -1347,6 +1400,8 @@ class SnapshotProcessor(HDF5Container):
for key in header: for key in header:
if units[header.index(key)] == "m": if units[header.index(key)] == "m":
df[key] = df[key] * self.info["unit_mass"].express(U.Msun) df[key] = df[key] * self.info["unit_mass"].express(U.Msun)
if "M" in df:
df["msink"] = df["M"]
else: else:
header = [ header = [
"id", "id",
@@ -1369,6 +1424,7 @@ class SnapshotProcessor(HDF5Container):
"Teff", "Teff",
] ]
df = pd.read_csv(csv_name, header=None, names=header) df = pd.read_csv(csv_name, header=None, names=header)
return {key: df[key].values for key in df} return {key: df[key].values for key in df}
def _pspec(self, overwrite_file=False, **kwargs): def _pspec(self, overwrite_file=False, **kwargs):
@@ -1566,7 +1622,7 @@ class SnapshotProcessor(HDF5Container):
# sort it and apply the sorting to the coordinates # 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... # this means that the particules of group 1 are written first then of group 2 etc...
ind_sort = np.argsort(group_ids) ind_sort = np.argsort(group_ids)
cells_group = val_mat[ind_sort] cells_group = cells_group[ind_sort]
cells_group[6] = group_ids[ind_sort] cells_group[6] = group_ids[ind_sort]
# Make it a pandas' DataFrame # Make it a pandas' DataFrame
@@ -1587,14 +1643,12 @@ class SnapshotProcessor(HDF5Container):
self.rules = { self.rules = {
# Base rules # Base rules
"coldens": Rule( "coldens": Rule(
self,
self._coldens, self._coldens,
"Column density", "Column density",
"/maps", "/maps",
unit=self.info["unit_density"] * self.info["unit_length"], unit=self.info["unit_density"] * self.info["unit_length"],
), ),
"T_mwavg": Rule( "T_mwavg": Rule(
self,
partial( partial(
self.ax_avg, self.ax_avg,
getter_T, getter_T,
@@ -1606,7 +1660,6 @@ class SnapshotProcessor(HDF5Container):
unit=self.info["unit_pressure"] / self.info["unit_density"], unit=self.info["unit_pressure"] / self.info["unit_density"],
), ),
"omega_mwavg": Rule( "omega_mwavg": Rule(
self,
partial(self._omega_average), partial(self._omega_average),
"Ax mass-weighted averaged azimuthal of angular frequency", "Ax mass-weighted averaged azimuthal of angular frequency",
"/maps", "/maps",
@@ -1614,7 +1667,6 @@ class SnapshotProcessor(HDF5Container):
/ (self.info["unit_length"] / self.lbox), / (self.info["unit_length"] / self.lbox),
), ),
"alpha_disk": Rule( "alpha_disk": Rule(
self,
self._alpha_disk, self._alpha_disk,
"Map of the Shakura&Sunaev alpha parameter for disks", "Map of the Shakura&Sunaev alpha parameter for disks",
"/maps", "/maps",
@@ -1627,7 +1679,6 @@ class SnapshotProcessor(HDF5Container):
], ],
), ),
"alpha_grav": Rule( "alpha_grav": Rule(
self,
self._alpha_grav, self._alpha_grav,
"Map of the graviational contrib to\ "Map of the graviational contrib to\
Shakura&Sunaev alpha parameter for disks", Shakura&Sunaev alpha parameter for disks",
@@ -1636,39 +1687,33 @@ class SnapshotProcessor(HDF5Container):
dependencies=["avg_map_coldens", "avg_map_T_mwavg"], dependencies=["avg_map_coldens", "avg_map_T_mwavg"],
), ),
"T": Rule( "T": Rule(
self,
self._temperature, self._temperature,
"Temperature slice", "Temperature slice",
"/maps", "/maps",
dependencies=["slice_rho"], dependencies=["slice_rho"],
unit=self.info["unit_pressure"] / self.info["unit_density"], unit=self.info["unit_pressure"] / self.info["unit_density"],
), ),
"levels": Rule( "levels": Rule(self._levels, "Max level within line of sight", "/maps"
self, self._levels, "Max level within line of sight", "/maps"
), ),
"jeans": Rule( "jeans": Rule(
self,
self._jeans, self._jeans,
"Jeans length slice", "Jeans length slice",
"/maps", "/maps",
dependencies=["slice_rho", "T"], dependencies=["slice_rho", "T"],
), ),
"jeans_ratio": Rule( "jeans_ratio": Rule(
self,
self._jeans_ratio, self._jeans_ratio,
"Jeans' length divided by the max resolution", "Jeans' length divided by the max resolution",
"/maps", "/maps",
dependencies=["jeans", "levels"], dependencies=["jeans", "levels"],
), ),
"Q": Rule( "Q": Rule(
self,
self._toomreQ_disk, self._toomreQ_disk,
"Toomre Q parameter for a Keplerian disk", "Toomre Q parameter for a Keplerian disk",
"/maps", "/maps",
dependencies=["coldens", "T_mwavg", "omega_mwavg"], dependencies=["coldens", "T_mwavg", "omega_mwavg"],
), ),
"sinks": Rule( "sinks": Rule(
self,
self._sinks, self._sinks,
group="/datasets", group="/datasets",
unit={ unit={
@@ -1702,20 +1747,17 @@ class SnapshotProcessor(HDF5Container):
"vz_gas": "unit_velocity", "vz_gas": "unit_velocity",
}, },
), ),
"pspec": Rule(self, self._pspec, "Power spectrum", "/hdf5"), "pspec": Rule(self._pspec, "Power spectrum", "/hdf5"),
"write_particles": Rule( "write_particles": Rule(self._write_particles, "Particles file", "/hdf5"
self, self._write_particles, "Particles file", "/hdf5"
), ),
"write_cells": Rule(self, self._write_cells, "Cells file", "/hdf5"), "write_cells": Rule(self._write_cells, "Cells file", "/hdf5"),
"filaments": Rule( "filaments": Rule(
self,
self._filaments, self._filaments,
"Filaments", "Filaments",
"/datasets", "/datasets",
dependencies={self.params.filaments.datamap: "z"}, dependencies={self.params.filaments.datamap: "z"},
), ),
"filaments_forces": Rule( "filaments_forces": Rule(
self,
self._filaments_forces, self._filaments_forces,
"Filaments", "Filaments",
"/datasets", "/datasets",
@@ -1730,14 +1772,12 @@ class SnapshotProcessor(HDF5Container):
), ),
# Helpers # Helpers
"radial_bins": Rule( "radial_bins": Rule(
self,
self._radial_bins, self._radial_bins,
"Radial bins", "Radial bins",
"/radial", "/radial",
unit=self.info["unit_length"] / self.lbox, unit=self.info["unit_length"] / self.lbox,
), ),
"radial_centers": Rule( "radial_centers": Rule(
self,
self._radial_centers, self._radial_centers,
"Centers of radial bins", "Centers of radial bins",
"/radial", "/radial",
@@ -1745,21 +1785,18 @@ class SnapshotProcessor(HDF5Container):
unit=self.info["unit_length"] / self.lbox, unit=self.info["unit_length"] / self.lbox,
), ),
"rr": Rule( "rr": Rule(
self,
self._rr, self._rr,
"Coordinate map", "Coordinate map",
"/maps", "/maps",
unit=self.info["unit_length"] / self.lbox, unit=self.info["unit_length"] / self.lbox,
), ),
"bins_on_map": Rule( "bins_on_map": Rule(
self,
self._bins_on_map, self._bins_on_map,
"Convert map coordinates to bins", "Convert map coordinates to bins",
"/maps", "/maps",
dependencies=["radial_bins", "rr"], dependencies=["radial_bins", "rr"],
), ),
"B_int": Rule( "B_int": Rule(
self,
self._B_int, self._B_int,
"Magnetic intensity slice", "Magnetic intensity slice",
"/maps", "/maps",
@@ -1768,14 +1805,12 @@ class SnapshotProcessor(HDF5Container):
), ),
# PDF # PDF
"rho_pdf": Rule( "rho_pdf": Rule(
self,
partial(self.vol_pdf, partial(simple_getter, "rho"), logbins=True), partial(self.vol_pdf, partial(simple_getter, "rho"), logbins=True),
"Global rho-PDF", "Global rho-PDF",
"/hist", "/hist",
unit=self.info["unit_density"], unit=self.info["unit_density"],
), ),
"rho_pdf_mw": Rule( "rho_pdf_mw": Rule(
self,
partial( partial(
self.vol_pdf, self.vol_pdf,
partial(simple_getter, "rho"), partial(simple_getter, "rho"),
@@ -1787,14 +1822,12 @@ class SnapshotProcessor(HDF5Container):
unit=self.info["unit_density"], unit=self.info["unit_density"],
), ),
"T_pdf": Rule( "T_pdf": Rule(
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"],
), ),
"cos_pdf": Rule( "cos_pdf": Rule(
self,
partial(self.cos_vfluct_B), partial(self.cos_vfluct_B),
"Global cos fluctuation-PDF", "Global cos fluctuation-PDF",
"/hist", "/hist",
@@ -1802,14 +1835,12 @@ class SnapshotProcessor(HDF5Container):
unit=U.none, unit=U.none,
), ),
"Brho": Rule( "Brho": Rule(
self,
self._Brho, self._Brho,
"Average of B as a function of rho", "Average of B as a function of rho",
"/datasets", "/datasets",
unit={"rho": self.info["unit_density"], "B": self.info["unit_mag"]}, unit={"rho": self.info["unit_density"], "B": self.info["unit_mag"]},
), ),
"Ek_Eb_rho": Rule( "Ek_Eb_rho": Rule(
self,
self._Ek_Eb_rho, self._Ek_Eb_rho,
"Average of Ek/Eb as a function of rho", "Average of Ek/Eb as a function of rho",
"/datasets", "/datasets",
@@ -1818,14 +1849,12 @@ class SnapshotProcessor(HDF5Container):
), ),
# Profiles # Profiles
"axis": Rule( "axis": Rule(
self,
partial(self._get_axis), partial(self._get_axis),
"Axis", "Axis",
"/profile", "/profile",
unit=self.info["unit_length"], unit=self.info["unit_length"],
), ),
"rho_prof": Rule( "rho_prof": Rule(
self,
partial(self._plane_avg_uniform, partial(simple_getter, "rho")), partial(self._plane_avg_uniform, partial(simple_getter, "rho")),
"Rho profile", "Rho profile",
"/profile", "/profile",
@@ -1834,28 +1863,24 @@ class SnapshotProcessor(HDF5Container):
), ),
# globals # globals
"time_num": Rule( "time_num": Rule(
self,
lambda: self.info["time"], lambda: self.info["time"],
"Time", "Time",
"/globals", "/globals",
unit=self.info["unit_time"], unit=self.info["unit_time"],
), ),
"mass": Rule( "mass": Rule(
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,
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"],
), ),
"mwa_sigma": Rule( "mwa_sigma": Rule(
self,
self._mwa_sigma, self._mwa_sigma,
"Mass weighted speed average", "Mass weighted speed average",
"/globals", "/globals",
@@ -1863,7 +1888,6 @@ class SnapshotProcessor(HDF5Container):
unit=self.info["unit_velocity"], unit=self.info["unit_velocity"],
), ),
"mwa_B_int": Rule( "mwa_B_int": Rule(
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",
@@ -1888,7 +1912,6 @@ class SnapshotProcessor(HDF5Container):
oct_getter = getter oct_getter = getter
self.rules["slice_" + name] = Rule( self.rules["slice_" + name] = Rule(
self,
partial(self.slice, getter, z=0.0, unit=unit), partial(self.slice, getter, z=0.0, unit=unit),
"{} slice".format(name), "{} slice".format(name),
"/maps", "/maps",
@@ -1896,7 +1919,6 @@ class SnapshotProcessor(HDF5Container):
) )
self.rules[name + "_mwavg"] = Rule( self.rules[name + "_mwavg"] = Rule(
self,
partial(self.ax_avg, oct_getter, mass_weighted=True, unit=unit), partial(self.ax_avg, oct_getter, mass_weighted=True, unit=unit),
"Ax mass-weighted averaged {}".format(name), "Ax mass-weighted averaged {}".format(name),
"/maps", "/maps",
@@ -1904,7 +1926,6 @@ class SnapshotProcessor(HDF5Container):
) )
self.rules[name + "_avg"] = Rule( self.rules[name + "_avg"] = Rule(
self,
partial(self.ax_avg, oct_getter, mass_weighted=False, unit=unit), partial(self.ax_avg, oct_getter, mass_weighted=False, unit=unit),
"Ax averaged {}".format(name), "Ax averaged {}".format(name),
"/maps", "/maps",
@@ -1944,7 +1965,8 @@ class SnapshotProcessor(HDF5Container):
# Norm # Norm
generic_rule( generic_rule(
field + "_norm", partial(norm_getter, field), self.unit_key[field] field + "_norm", partial(norm_getter, field), self.unit_key[field],
oct_getter=partial(oct_norm_getter, field)
) )
else: else:
@@ -1958,7 +1980,6 @@ class SnapshotProcessor(HDF5Container):
self.rules["rad_avg_" + name] = Rule( self.rules["rad_avg_" + name] = Rule(
self,
partial(self._rad_avg, name), partial(self._rad_avg, name),
"Azimuthal average of {}".format(name), "Azimuthal average of {}".format(name),
"/radial", "/radial",
@@ -1966,7 +1987,6 @@ class SnapshotProcessor(HDF5Container):
unit=unit, unit=unit,
) )
self.rules["rad_mwavg_" + name] = Rule( self.rules["rad_mwavg_" + name] = Rule(
self,
partial(self._rad_avg, name, mass_weighted=True), partial(self._rad_avg, name, mass_weighted=True),
"Mass weighted azimuthal average of {}".format(name), "Mass weighted azimuthal average of {}".format(name),
"/radial", "/radial",
@@ -1974,7 +1994,6 @@ class SnapshotProcessor(HDF5Container):
unit=unit, unit=unit,
) )
self.rules["rad_avg_map_" + name] = Rule( self.rules["rad_avg_map_" + name] = Rule(
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),
"/maps", "/maps",
@@ -1982,7 +2001,6 @@ class SnapshotProcessor(HDF5Container):
unit=unit, unit=unit,
) )
self.rules["rad_mwavg_map_" + name] = Rule( self.rules["rad_mwavg_map_" + name] = Rule(
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),
"/maps", "/maps",
@@ -1990,14 +2008,12 @@ class SnapshotProcessor(HDF5Container):
unit=unit, unit=unit,
) )
self.rules["fluct_" + name] = Rule( self.rules["fluct_" + name] = Rule(
self,
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, "rad_avg_map_" + name], dependencies=[name, "rad_avg_map_" + name],
) )
self.rules["dispersion_rad_" + name] = Rule( self.rules["dispersion_rad_" + name] = Rule(
self,
partial(self._dispersion_rad, name), partial(self._dispersion_rad, name),
"radial RMS of {}".format(name), "radial RMS of {}".format(name),
"/radial", "/radial",
@@ -2005,21 +2021,18 @@ class SnapshotProcessor(HDF5Container):
unit=unit, unit=unit,
) )
self.rules["mwfluct_" + name] = Rule( self.rules["mwfluct_" + name] = Rule(
self,
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, "rad_mwavg_map_" + name], dependencies=[name, "rad_mwavg_map_" + name],
) )
self.rules["fluct_pdf_" + name] = Rule( self.rules["fluct_pdf_" + name] = Rule(
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["fluct_mwpdf_" + name] = Rule( self.rules["fluct_mwpdf_" + name] = Rule(
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(
name name
@@ -2028,7 +2041,6 @@ class SnapshotProcessor(HDF5Container):
dependencies=["rr", "mwfluct_" + name], dependencies=["rr", "mwfluct_" + name],
) )
self.rules["fit_fluct_pdf_" + name] = Rule( self.rules["fit_fluct_pdf_" + name] = Rule(
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",
@@ -2038,7 +2050,6 @@ class SnapshotProcessor(HDF5Container):
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,
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",
studyprocessor.py
+32 -25
View File
@@ -63,10 +63,13 @@ class StudyProcessor(Aggregator, HDF5Container):
self.runs = selector.runs self.runs = selector.runs
self.nums = selector.nums self.nums = selector.nums
run0 = self.runs[0]
self.info = selector.info[run0][self.nums[run0][0]]
self.namelist = selector.namelist
# Get postprocesor objets for each run and infos on them # Get postprocesor objets for each run and infos on them
self.snaps = {} self.snaps = {}
self.info = {}
for run in self.runs: for run in self.runs:
path_run = path + "/" + run path_run = path + "/" + run
path_out_run = self.path_out + "/" + run path_out_run = self.path_out + "/" + run
@@ -82,9 +85,7 @@ class StudyProcessor(Aggregator, HDF5Container):
unit_time=unit_time, unit_time=unit_time,
) )
run0 = self.runs[0]
self.info = selector.info[run0][self.nums[run0][0]]
self.namelist = selector.namelist
# Save namelist and logs # Save namelist and logs
if self.params.out.copy_info: if self.params.out.copy_info:
@@ -234,8 +235,19 @@ class StudyProcessor(Aggregator, HDF5Container):
value = snap.get_value(node_name) value = snap.get_value(node_name)
return value return value
def get_nml(self, nml_key, run): def get_nml(self, nml_key=None, run=None):
return self.namelist[run][nml_key] if run is not None:
if nml_key is not None:
return self.namelist[run][nml_key]
else:
return self.namelist[run]
else:
if nml_key is not None:
return {run : self.namelist[run][nml_key] for run in self.runs}
else:
return self.namelist
def get_pdf_slope(self, name, run, num): def get_pdf_slope(self, name, run, num):
snap = self.snaps[run][num] snap = self.snaps[run][num]
@@ -464,15 +476,19 @@ class StudyProcessor(Aggregator, HDF5Container):
time_gas = self.get_value("/series/coarse_step_from_log/time") time_gas = self.get_value("/series/coarse_step_from_log/time")
mass_gas = self.get_value("/series/coarse_step_from_log/mcons") mass_gas = self.get_value("/series/coarse_step_from_log/mcons")
mass_sink = self.get_value("/series/sinks_from_log/mass_sink") mass_sink = self.get_value("/series/sinks_from_log/mass_sink")
time_sink = self.get_value("/series/sinks_from_log/time")
total_mass = dict() total_mass = dict()
for run in self.runs: for run in self.runs:
if time_sink[run][-1] > time_gas[run][-1]:
time_sink[run] = time_sink[run][:-1]
mass_sink[run] = mass_sink[run][:-1]
# A bit specific ... needs to be generalized (TODO) # A bit specific ... needs to be generalized (TODO)
info = self.snaps[run][self.nums[run][0]].info info = self.snaps[run][self.nums[run][0]].info
surface = (info["unit_length"].express(U.pc)) ** 2 surface = (info["unit_length"].express(U.pc)) ** 2
m0 = self.get_coldens0(run) * surface # Initial mass in Msun m0 = self.get_coldens0(run) * surface # Initial mass in Msun
offset = time_gas[run].size - mass_sink[run].size offset = time_gas[run].size - time_sink[run].size
mass_gas[run] = m0 + m0*mass_gas[run] # convert in Msun mass_gas[run] = m0 + m0*mass_gas[run] # convert in Msun
total_mass[run] = mass_gas[run].copy() total_mass[run] = mass_gas[run].copy()
total_mass[run][offset:] = mass_gas[run][offset:] + mass_sink[run] # re add sink_mass total_mass[run][offset:] = mass_gas[run][offset:] + mass_sink[run] # re add sink_mass
@@ -574,7 +590,6 @@ class StudyProcessor(Aggregator, HDF5Container):
name = "time_" + glob_name name = "time_" + glob_name
self.rules[name] = Rule( self.rules[name] = Rule(
self,
partial( partial(
self._time_series, self._time_series,
partial( partial(
@@ -627,7 +642,6 @@ class StudyProcessor(Aggregator, HDF5Container):
) )
self.rules[name] = Rule( self.rules[name] = Rule(
self,
fn, fn,
group="/comp", group="/comp",
unit=units, unit=units,
@@ -640,7 +654,6 @@ class StudyProcessor(Aggregator, HDF5Container):
self.rules = { self.rules = {
# Read from log # Read from log
"sinks_from_log": Rule( "sinks_from_log": Rule(
self,
partial( partial(
self._from_log, self._from_log,
["time", "mass_sink", "nb_sink"], ["time", "mass_sink", "nb_sink"],
@@ -655,7 +668,6 @@ class StudyProcessor(Aggregator, HDF5Container):
}, },
), ),
"issfr": Rule( "issfr": Rule(
self,
self._ssfr_from_mass_sink, self._ssfr_from_mass_sink,
group="/series/sinks_from_log", group="/series/sinks_from_log",
unit=U.ssfr, unit=U.ssfr,
@@ -663,7 +675,6 @@ class StudyProcessor(Aggregator, HDF5Container):
dependencies=["sinks_from_log"], dependencies=["sinks_from_log"],
), ),
"ssm": Rule( "ssm": Rule(
self,
self._surfacic_sink_mass, self._surfacic_sink_mass,
group="/series/sinks_from_log", group="/series/sinks_from_log",
unit=U.Msun / U.pc ** 2, unit=U.Msun / U.pc ** 2,
@@ -671,7 +682,6 @@ class StudyProcessor(Aggregator, HDF5Container):
dependencies=["sinks_from_log"], dependencies=["sinks_from_log"],
), ),
"sfr_from_log": Rule( "sfr_from_log": Rule(
self,
partial(self._from_log, ["time", "sfr"], self._extract_sfr_from_log), partial(self._from_log, ["time", "sfr"], self._extract_sfr_from_log),
group="/series", group="/series",
unit={"time": U.year, "sfr": U.ssfr}, unit={"time": U.year, "sfr": U.ssfr},
@@ -681,7 +691,6 @@ class StudyProcessor(Aggregator, HDF5Container):
}, },
), ),
"rms_from_log": Rule( "rms_from_log": Rule(
self,
partial( partial(
self._from_log, self._from_log,
["time", "dt", "turb_rms", "turb_energy"], ["time", "dt", "turb_rms", "turb_energy"],
@@ -706,7 +715,6 @@ class StudyProcessor(Aggregator, HDF5Container):
}, },
), ),
"coarse_step_from_log": Rule( "coarse_step_from_log": Rule(
self,
partial( partial(
self._from_log, self._from_log,
[ [
@@ -738,7 +746,6 @@ class StudyProcessor(Aggregator, HDF5Container):
}, },
), ),
"fine_step_from_log": Rule( "fine_step_from_log": Rule(
self,
partial( partial(
self._from_log, self._from_log,
["time", "fine_step", "dt", "a", "mem_cells", "mem_parts"], ["time", "fine_step", "dt", "a", "mem_cells", "mem_parts"],
@@ -755,7 +762,6 @@ class StudyProcessor(Aggregator, HDF5Container):
}, },
), ),
"stellar_from_log": Rule( "stellar_from_log": Rule(
self,
partial( partial(
self._from_log, self._from_log,
["time", "mass", "lifetime", "id"], ["time", "mass", "lifetime", "id"],
@@ -770,7 +776,6 @@ class StudyProcessor(Aggregator, HDF5Container):
}, },
), ),
"turb_power": Rule( "turb_power": Rule(
self,
self._turb_power, self._turb_power,
group="/series/rms_from_log", group="/series/rms_from_log",
unit={ unit={
@@ -784,7 +789,6 @@ class StudyProcessor(Aggregator, HDF5Container):
), ),
# Read from outputs # Read from outputs
"time": Rule( "time": Rule(
self,
partial( partial(
self._time_series, partial(self.get_global, "/globals/time_num") self._time_series, partial(self.get_global, "/globals/time_num")
), ),
@@ -793,7 +797,6 @@ class StudyProcessor(Aggregator, HDF5Container):
dependencies=["time_num"], dependencies=["time_num"],
), ),
"time_rho_prof": Rule( "time_rho_prof": Rule(
self,
partial( partial(
self._time_series, partial(self.get_snap_value, "/profile/rho_prof") self._time_series, partial(self.get_snap_value, "/profile/rho_prof")
), ),
@@ -801,15 +804,20 @@ class StudyProcessor(Aggregator, HDF5Container):
dependencies={"time": None, "rho_prof": "__parent__"}, dependencies={"time": None, "rho_prof": "__parent__"},
), ),
"time_coldens_pdf": Rule( "time_coldens_pdf": Rule(
self,
partial( partial(
self._time_series, partial(self.get_snap_value, "/hist/pdf_coldens") self._time_series, partial(self.get_snap_value, "/hist/pdf_coldens")
), ),
group="/series", group="/series",
dependencies={"time": None, "pdf_coldens": "__parent__"}, dependencies={"time": None, "pdf_coldens": "__parent__"},
), ),
"time_rho_pdf": Rule(
partial(
self._time_series, partial(self.get_snap_value, "/hist/rho_pdf")
),
group="/series",
dependencies={"time": None},
),
"time_pdf_slope_coldens": Rule( "time_pdf_slope_coldens": Rule(
self,
partial( partial(
self._time_series, self._time_series,
partial( partial(
@@ -822,7 +830,6 @@ class StudyProcessor(Aggregator, HDF5Container):
dependencies={"time": None, "fit_pdf_coldens": "z"}, dependencies={"time": None, "fit_pdf_coldens": "z"},
), ),
"sbeta_onavg": Rule( "sbeta_onavg": Rule(
self,
partial(self._sbeta_onavg), partial(self._sbeta_onavg),
group="/comp", group="/comp",
dependencies={ dependencies={
@@ -832,7 +839,6 @@ class StudyProcessor(Aggregator, HDF5Container):
), ),
# namelist # namelist
"nml": Rule( "nml": Rule(
self,
lambda nml_key: self._compare( lambda nml_key: self._compare(
partial(self.get_nml, nml_key), use_num=False partial(self.get_nml, nml_key), use_num=False
), ),
@@ -854,6 +860,7 @@ class StudyProcessor(Aggregator, HDF5Container):
"turb_power", "turb_power",
"time_rho_prof", "time_rho_prof",
"time_coldens_pdf", "time_coldens_pdf",
"time_rho_pdf"
]: ]:
self._gen_rule_avg(name) self._gen_rule_avg(name)