diff --git a/__init__.py b/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/__pycache__/comparator.cpython-37.pyc b/__pycache__/comparator.cpython-37.pyc deleted file mode 100644 index e38ea18..0000000 Binary files a/__pycache__/comparator.cpython-37.pyc and /dev/null differ diff --git a/aggregator.py b/aggregator.py index 1d43d2a..e3f349d 100644 --- a/aggregator.py +++ b/aggregator.py @@ -4,7 +4,7 @@ import numpy as np from functools import partial from baseprocessor import Rule import snapshotprocessor -from mypool import MyPool +from utils.mypool import MyPool try: from mpi4py.futures import MPIPoolExecutor diff --git a/baseprocessor.py b/baseprocessor.py index 1fb8e9b..1f21f50 100644 --- a/baseprocessor.py +++ b/baseprocessor.py @@ -10,11 +10,11 @@ from functools import partial import numpy as np import tables -from tables import HDF5ExtError +from tables import HDF5ExtError, NoSuchNodeError from tables.registry import class_name_dict -from params import default_params, load_params -from units import U +from utils.params import default_params, load_params +from utils.units import U import sys import traceback @@ -44,6 +44,7 @@ class Rule: return self.process_fn(arg, **kwargs) else: return self.process_fn(**kwargs) + class BaseProcessor: """ Base class for processors, should not be instanciated @@ -292,8 +293,17 @@ class HDF5Container(BaseProcessor): ) else: value = node.read() + if isinstance(unit, dict): + name = os.path.basename(node_name) + if name in unit: + unit = unit[name] + else: + unit = None if not (unit is None or unit_old is None or unit_old == U.none): value = value * unit_old.express(unit) + except NoSuchNodeError: + self.logger.error(f"The value {node_name} is node available", stack_info=True) + raise finally: if not open_before: self.close() diff --git a/fragdisk/__init__.py b/fragdisk/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/gui.py b/fragdisk/gui.py similarity index 99% rename from gui.py rename to fragdisk/gui.py index 6561773..52055ec 100644 --- a/gui.py +++ b/fragdisk/gui.py @@ -17,7 +17,7 @@ from scipy.stats import linregress from skimage.draw import line from snapshotprocessor import SnapshotProcessor -from params import default_params +from utils.params import default_params class DraggablePoint: diff --git a/galactica/__init__.py b/galactica/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/fragdisk_galactica.py b/galactica/fragdisk_galactica.py similarity index 100% rename from fragdisk_galactica.py rename to galactica/fragdisk_galactica.py diff --git a/ismfeed_galactica.py b/galactica/ismfeed_galactica.py similarity index 100% rename from ismfeed_galactica.py rename to galactica/ismfeed_galactica.py diff --git a/ramses_astrophysix.py b/galactica/ramses_astrophysix.py similarity index 100% rename from ramses_astrophysix.py rename to galactica/ramses_astrophysix.py diff --git a/galturb/__init__.py b/galturb/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/galaxy.py b/galturb/galaxy.py similarity index 100% rename from galaxy.py rename to galturb/galaxy.py diff --git a/sectors_extraction.py b/galturb/sectors_extraction.py similarity index 100% rename from sectors_extraction.py rename to galturb/sectors_extraction.py diff --git a/ism.py b/ism.py deleted file mode 100644 index 2fa0ddf..0000000 --- a/ism.py +++ /dev/null @@ -1,101 +0,0 @@ -# coding: utf-8 - -import numpy as np -import pandas as pd -from plotter import U -import snapshotprocessor - -mp = 1.4 * 1.66 * 10 ** (-24) * U.g -z0 = 150 * U.pc -sink_header = [ - "Id", - "M", - "dmf", - "x", - "y", - "z", - "vx", - "vy", - "vz", - "rot_period", - "lx", - "ly", - "lz", - "acc_rate", - "acc_lum", - "age", - "int_lum", - "Teff", -] - - -def convert_coldens_s(n0): - return (np.sqrt(2 * np.pi) * mp * z0 * (n0 * U.cm ** (-3))).express(U.coldens) - - -convert_coldens = np.vectorize(convert_coldens_s) - - -def get_dispersion(dset, name): - """ - Compute dispersion from dset[name] - """ - vel = dset[name] - mass = dset["mass"] - mass_tot = np.sum(mass) - if mass_tot > 0: - mean = np.sum(mass * vel) / mass_tot - return np.sqrt(np.sum(mass * (vel - mean) ** 2) / mass_tot) - else: - return 0 - -def get_sinks(pp): - csv_name = f"{pp.path}/output_{pp.num:05}/sink_{pp.num:05}.csv" - return pd.read_csv(csv_name, header=None, names=sink_header) - - -def analyze_box(pp): - pp.cells["mass"] = snapshotprocessor.mass_func(pp.cells) - coldens = pp.get_value("/maps/coldens_z", unit=U.coldens) - sinks = get_sinks(pp) - sinks["mass"] = sinks.M - res = {} - dirs = ["x", "y", "z"] - res["time"] = pp.info["time"] * pp.info["unit_time"].express(U.Myr) - for i, dir in enumerate(dirs): - pp.cells[f"v{dir}"] = pp.cells["vel"][:, i] - res[f"sigma_{dir}"] = get_dispersion(pp.cells, f"v{dir}") * pp.info[ - "unit_velocity" - ].express(U.km_s) - res[f"sigma_sinks_{dir}"] = get_dispersion(sinks, f"v{dir}") * pp.info[ - "unit_velocity" - ].express(U.km_s) - res["coldens_mean"] = np.mean(coldens) - res["n0"] = pp.get_nml("cloud_params/dens0") - res["mass"] = np.sum( - pp.cells["mass"] - * (pp.info["unit_density"] * pp.info["unit_length"] ** 3).express(U.Msun) - ) - res["coldens_initial"] = convert_coldens_s(res["n0"]) - res["mass_initial"] = res["coldens_initial"] * 1e6 - res["coldens_mean"] = np.mean(coldens) - res["coldens_beam"] = res["mass"] / (pp.info["unit_length"].express(U.pc)) ** 2 - res["nsink"] = sinks.M.count() - res["mass_sink"] = np.sum(sinks.M) - - return res - - -def load_wrapper(pp, fun): - """ - Wrapper to load_unload data and map function - """ - pp.load_cells(keys=["pos", "vel", "dx", "rho"]) - pp.coldens("z") - res = fun(pp) - pp.unload_cells() - return res - - -def allinone(pp): - return load_wrapper(pp, analyze_box) diff --git a/ismfeed/__init__.py b/ismfeed/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/ismfeed.py b/ismfeed/ismfeed.py similarity index 100% rename from ismfeed.py rename to ismfeed/ismfeed.py diff --git a/pipeline.py b/pipeline.py index 54d608d..e6e8871 100644 --- a/pipeline.py +++ b/pipeline.py @@ -9,7 +9,7 @@ import numpy as np from plotter import Plotter, plt from snapshotprocessor import SnapshotProcessor, get_time from studyprocessor import StudyProcessor -from params import default_params, load_params +from utils.params import default_params, load_params fake_pp = SnapshotProcessor() diff --git a/plotter.py b/plotter.py index 41eda37..0b387fd 100644 --- a/plotter.py +++ b/plotter.py @@ -36,12 +36,11 @@ except ModuleNotFoundError: print("WARNING: no movie support (missing module moviepy)") from mpl_toolkits.axes_grid1.inset_locator import inset_axes -import pspec_read from baseprocessor import Rule, BaseProcessor from aggregator import Aggregator from studyprocessor import StudyProcessor -from run_selector import RunSelector -from units import U, unit_str, convert_exp +from runselector import RunSelector +from utils.units import U, unit_str, convert_exp try: import lic @@ -56,7 +55,7 @@ from astrophysix.simdm.experiment import ( ParameterVisibility, Simulation, ) -from ramses_astrophysix import ramses +from galactica.ramses_astrophysix import ramses filetype_from_ext = {ext: ft for ft in FileType for ext in ft.extension_list} @@ -270,6 +269,7 @@ class Plotter(Aggregator, BaseProcessor): self.nums, self.path_out, self.params, + tag=tag, unit_time=unit_time, selector=self.selector, ) @@ -1434,15 +1434,6 @@ class Plotter(Aggregator, BaseProcessor): if subname_y: hdf5_y.close() - def _pspec(self, name, **kwargs): - """ - Plot power spectrum (wrapper around pspec_read) - """ - del kwargs["run"] - file_pspec = self.current_processor.get_value("/hdf5/pspec") - num = self.current_processor.num - getattr(pspec_read, "pspec_" + name)(file_pspec, ".", num, **kwargs) - def _overlay_fit(self, x, y, yerr=None, kind="linear", label=None, **kwargs): """ Add an overlay : fit a curve, linear or powerlaw @@ -1527,7 +1518,7 @@ class Plotter(Aggregator, BaseProcessor): This is where rules are defined """ self.rules = { - "plot_arrays": PlotRule(lambda arg, **kwargs: self._plot(*arg, **kwargs), kind="comp" + "plot_comp": PlotRule(lambda arg, **kwargs: self._plot(*arg, **kwargs), kind="comp" ), "plot_run": PlotRule(lambda arg, **kwargs: self._plot(*arg, **kwargs), kind="run" ), @@ -1674,7 +1665,6 @@ class Plotter(Aggregator, BaseProcessor): "Density profile", dependencies=["axis", "rho_prof"], ), - "pspec": PlotRule(self, self._pspec, dependencies={"pspec": None}), "sbeta": PlotRule( partial( self._plot, @@ -1965,13 +1955,15 @@ class Plotter(Aggregator, BaseProcessor): self._gen_from_log("fine_step_from_log", name) for name in ["time", "dt", "a", "mem_cells", "mem_parts"]: self._gen_from_log("fine_step_from_log", name_y=name, name_x="fine_step") + + self._gen_from_log("SN_momentum_from_log", name_x="time", name_y="SN_momentum") # Dict of overlays self.overlays = { "g": partial(self._overlay_vector, "g"), "B": self._overlay_B, "vel": self._overlay_speed, - "speed": self._overlay_speed, + "speed": self._overlay_speed, "levels": self._overlay_levels, "contour": self._overlay_contour, "particles": self._overlay_particles, diff --git a/pspec_read.py b/pspec_read.py deleted file mode 100644 index 7f17211..0000000 --- a/pspec_read.py +++ /dev/null @@ -1,374 +0,0 @@ -from builtins import str - -import matplotlib.pyplot as P -import numpy as np -import tables as T - -################################################################ -# 3D -################################################################ - - -def pspec_rho(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/rho") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 2) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^2 \mathcal{P}(\rho)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_rho_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_rho_" + format(num, "05") + ".jpeg") - - -def pspec_B(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/B") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 2 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^2 \mathcal{P}(B)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_B_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_B_" + format(num, "05") + ".jpeg") - - -def pspec_v(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/v") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 4 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^4 \mathcal{P}(v)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_v_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_v_" + format(num, "05") + ".jpeg") - - -def pspec_vz(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/vz") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 4 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^4 \mathcal{P}(vz)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_vz_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_vz_" + format(num, "05") + ".jpeg") - - -def pspec_cos(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/cos_vB") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 2 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^2 \mathcal{P}(\cos(\vec{v}.\vec{B}))$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_cos_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_cos_" + format(num, "05") + ".jpeg") - - -def pspec_logrho(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/logrho") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 2 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^2 \mathcal{P}(\log(\rho))$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_logrho_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_logrho_" + format(num, "05") + ".jpeg") - - -def pspec_kr(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/kr") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 4) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^4 \mathcal{P}(\rho^{1/3}v)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_kr_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_kr_" + format(num, "05") + ".jpeg") - - -def pspec_vc(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/vc") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 4) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^4 \mathcal{P}(v_c)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_vc_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_vc_" + format(num, "05") + ".jpeg") - - -def pspec_vs(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/vs") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 4) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^4 \mathcal{P}(v_s)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_vs_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_vs_" + format(num, "05") + ".jpeg") - - -################################################################### -# 2D -################################################################### - - -def pspec_rho_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/rho") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 1) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k \mathcal{P}(\rho)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_rho_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_rho_2D_" + format(num, "05") + ".jpeg") - - -def pspec_B_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/B") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 1 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k \mathcal{P}(B)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_B_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_B_2D_" + format(num, "05") + ".jpeg") - - -def pspec_v_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/v") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 3 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^3 \mathcal{P}(v)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_v_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_v_2D_" + format(num, "05") + ".jpeg") - - -def pspec_vz_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/vz") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 3 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^3 \mathcal{P}(v_z)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_vz_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_vz_2D_" + format(num, "05") + ".jpeg") - - -def pspec_cos_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/cos_vB") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k \mathcal{P}(\cos(\vec{v}.\vec{B}))$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_cos_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_cos_2D_" + format(num, "05") + ".jpeg") - - -def pspec_logrho_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/logrho") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k \mathcal{P}(\log(\rho))$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_logrho_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_logrho_2D_" + format(num, "05") + ".jpeg") - - -def pspec_kr_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/kr") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 3) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^3 \mathcal{P}(\rho^{1/3}v)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_kr_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_kr_2D_" + format(num, "05") + ".jpeg") - - -def pspec_vc_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/vc") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 3) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^3 \mathcal{P}(v_c)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_vc_2D" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_vc_2D" + format(num, "05") + ".jpeg") - - -def pspec_vs_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/vs") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 3) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^3 \mathcal{P}(v_s)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_vs_2D" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_vs_2D" + format(num, "05") + ".jpeg") diff --git a/run_selector.py b/runselector.py similarity index 100% rename from run_selector.py rename to runselector.py diff --git a/scripts/__init__.py b/scripts/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/pspec_new.py b/scripts/pspec.py similarity index 100% rename from pspec_new.py rename to scripts/pspec.py diff --git a/select_snapshot.py b/select_snapshot.py index 032a7f8..d954af3 100644 --- a/select_snapshot.py +++ b/select_snapshot.py @@ -5,7 +5,7 @@ Select snaphots with a criterion """ -from run_selector import RunSelector +from runselector import RunSelector from plotter import Plotter, U import os diff --git a/snapshotprocessor.py b/snapshotprocessor.py index fc245d9..e8e3474 100644 --- a/snapshotprocessor.py +++ b/snapshotprocessor.py @@ -34,10 +34,10 @@ from pymses.sources.hop.hop import HOP from fil_finder import FilFinder2D from scipy import fft -import pspec_new +from scripts import pspec as pspec -from units import U +from utils.units import U from baseprocessor import ( HDF5Container, Rule, @@ -48,7 +48,7 @@ from baseprocessor import ( oct_vect_getter, ) -from run_selector import RunSelector +from runselector import RunSelector # Getters @@ -212,10 +212,10 @@ def pspec(map2D): # Compute fft fmap = fft.fft2(map2D) # Compute the power map from the fft - pmap = pspec_new.pcube(fmap) + pmap = pspec.pcube(fmap) # Use the power map and the fft to compute the powerspectrum # This is typically an histogram of k weighted by the fourier transform value - pspec, kbins, pspec2, fbins = pspec_new.pspectrum(pmap, kmap, kbins, 1, 0) + pspec, kbins, pspec2, fbins = pspec.pspectrum(pmap, kmap, kbins, 1, 0) # Return bin center and power spectrum return 0.5 * (kbins[1:] + kbins[:-1]), pspec @@ -1434,10 +1434,10 @@ class SnapshotProcessor(HDF5Container): return {key: df[key].values for key in df} - def _pspec(self, overwrite_file=False, **kwargs): + def pspec(self, overwrite_file=False, **kwargs): outfile = self.pspec_filename if overwrite_file or not os.path.exists(self.pspec_filename): - pspec_new.pspec(repo=self.path, iouts=[self.num], outfile=outfile, **kwargs) + pspec.pspec(repo=self.path, iouts=[self.num], outfile=outfile, **kwargs) return np.array([self.pspec_filename]) def _write_particles(self): @@ -1754,7 +1754,7 @@ class SnapshotProcessor(HDF5Container): "vz_gas": "unit_velocity", }, ), - "pspec": Rule(self._pspec, "Power spectrum", "/hdf5"), + "pspec": Rule(self.pspec, "Power spectrum", "/hdf5"), "write_particles": Rule(self._write_particles, "Particles file", "/hdf5" ), "write_cells": Rule(self._write_cells, "Cells file", "/hdf5"), diff --git a/studyprocessor.py b/studyprocessor.py index d1f0c7b..ffa48a5 100644 --- a/studyprocessor.py +++ b/studyprocessor.py @@ -12,9 +12,9 @@ from scipy.stats import linregress from baseprocessor import Rule, HDF5Container from aggregator import Aggregator from snapshotprocessor import SnapshotProcessor -from run_selector import RunSelector -from params import default_params -from units import U +from runselector import RunSelector +from utils.params import default_params +from utils.units import U class StudyProcessor(Aggregator, HDF5Container): @@ -452,6 +452,31 @@ class StudyProcessor(Aggregator, HDF5Container): series["turb_energy"][run].append(np.nan) return series + def _extract_SN_Mom_from_log(self, series, log_filename, run): + cmd_grep = f"grep -e 'SN momentum' -e 'Fine step' {log_filename}" + content = os.popen(cmd_grep).readlines() + block_err = [] + time = 0 + for i in range(0, len(content)): + try: + if content[i][1:5] == "Fine": + data = content[i].replace("=", " ").split() + time = np.float(data[4]) + elif content[i][1:3] == "SN" : + series["time"][run].append(time) + series["SN_momentum"][run].append(np.float(content[i].split()[-1])) + else: + raise ValueError("Wrong start of line") + except (AssertionError, ValueError, IndexError): + block_err.append(i) + + if len(block_err) > 0: + self.logger.warning( + f"Error encountered in parsing {log_filename} (grepped blocks {block_err})" + ) + return series + + def get_logs(self, run): glob_str = f"{self.path}/{run}/{self.params.input.log_prefix}*" logs = glob.glob(glob_str) @@ -815,6 +840,15 @@ class StudyProcessor(Aggregator, HDF5Container): "id": U.none, }, ), + "SN_momentum_from_log": Rule( + partial(self._from_log, ["time", "SN_momentum"], self._extract_SN_Mom_from_log), + group="/series", + unit={"time": "unit_time", "SN_momentum" : {"unit_mass" : 1, "unit_velocity" : 1}}, + description={ + "time": "Time", + "SN_momentum": "Injected momentum", + }, + ), "turb_power": Rule( self._turb_power, group="/series/rms_from_log", diff --git a/turbox/__init__.py b/turbox/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/turbox/turbox.py b/turbox/turbox.py new file mode 100644 index 0000000..fdb8038 --- /dev/null +++ b/turbox/turbox.py @@ -0,0 +1,160 @@ +""" + Toolbox for the TURBOX project (simulation part). + + By Noé Brucy and Corentin Le Yuehlic, 2022 + + Compute and plot eta_d, from the slope of the logdensity power spectrum +""" + +import numpy as np +from plotter import U +from baseprocessor import Rule +from matplotlib import pyplot as plt +import pandas as pd +import tables +from scipy.stats import linregress + + +def get_pspec(pp, field:str, dim:int=3): + """Read power spectruù + + Parameters + ---------- + pp : SnapshotProcessor + field : str + field to get + dim : int, optional + dimension (2 or 3), by default 3 + + Returns + ------- + tupple (np.array, np.array) + wave number and correponding powers + """ + h5file = tables.File(pp.pspec_filename) + path = f"/out_{pp.num:05}/d{dim}/{field}" + node = h5file.get_node(path) + kbins = node.kbins.read() + pspec = node.pspec.read() + h5file.close() + k = (kbins[:-1] + kbins[1:]) / 2 + return k, pspec + + +span_resolution = { + 256: (0.8, 1.1), + 512: (0.5, 1.7), + 1024: (0.4, 1.7) +} + + +def get_slope(pp, field:str, resol:int, plotdebug:bool=False): + """Get the slope of the Power specturm using linear regression in the selected range + + Parameters + ---------- + pp : Snapshotprocessor + field : str + field name + resol : int + resolution (number of cells on 1 side of the simulation cube) + + Returns + ------- + tuple (float, float) + Slope, square value of the correlation coefficient + """ + # Trustworthy span od the power spectrum in log10(k) as a function of the resolution + + logkmin, logkmax = span_resolution[resol] + k, power = get_pspec(pp, field) + logk, logpower = np.log10(k), np.log10(power) + mask = (logk >= logkmin) & (logk < logkmax) + results = linregress(logk[mask], logpower[mask]) + if plotdebug: + plt.figure() + plt.plot(logk, logpower) + plt.plot(logk[mask], results.slope*logk[mask]+ results.intercept, lw=3, ls=":", color="k") + pp.logger.info(f"Fit results in get_slope({field}, {resol}): slope:{results.slope:.2f}, b:{results.intercept:.2f}, R2:{results.rvalue**2:.2f}") + if results.rvalue**2 < 0.8: + pp.logger.warning(f"Bad fit in get_slope({field}, {resol}) with {logkmin} <= logk < {logkmax}") + pp.logger.warning(f"log(k) is \n {logk[mask]}") + pp.logger.warning(f"log(power) is \n {logpower[mask]}") + + return results.slope, results.intercept, results.rvalue**2 + +def build_suite(pl, redo=False, cs0=0.28834810480560674): + """Compute an array of parameter for each run in the Plotter pl + + Parameters + ---------- + pl : Plotter + redo : bool, optional + redo the comptutation of the velocity dispersion, by default False + cs0 : float, optional + Sound speed in the simulations, by default 0.28834810480560674 + + Returns + ------- + pd.Dataframe + dataframe with the properties of the simulation + """ + + df = dict() + df["snapshots"] = pl.nums.values() + df["n0"] = pl.study.get_nml("galbox_params/dens0").values() + df["turbinit"] = pl.study.get_nml("galbox_params/turb").values() + df["solver"] = pl.study.get_nml("hydro_params/riemann").values() + df["slope"] = pl.study.get_nml("hydro_params/slope_type").values() + df["res"] = pl.study.get_nml("amr_params/levelmin").values() + df["frms"] = pl.study.get_nml("turb_params/turb_rms").values() + df["seed"] = pl.study.get_nml("turb_params/turb_seed").values() + + df["comp"] = pl.study.get_nml("turb_params/comp_frac").values() + df["L"] = pl.study.get_nml("amr_params/boxlen").values() + + df["T_turb"] = (np.array(list(pl.study.get_nml("turb_params/turb_T").values())) + * pl.study.info["unit_time"].express(U.Myr)) + df = pd.DataFrame(df, index=pl.runs) + + if redo: + pl.study.avg_time_sigma("x", overwrite_dep=False) + pl.study.avg_time_sigma("y", overwrite_dep=False) + pl.study.avg_time_sigma("z", overwrite_dep=False) + pl.study.time(overwrite=True) + + for ax in ["x", "y", "z"]: + df[f"sigma_{ax}"] = np.array(list(map( + lambda x : x.T[0], + [pl.study.get_value(f"/series/time_sigma_{ax}", + unit=U.km_s)[run] for run in pl.runs]))) + + df["sigma_all"] = df[f"sigma_x"]**2 + df[f"sigma_y"]**2 + df[f"sigma_z"]**2 + df["sigma_all"] = list(map(np.sqrt, df["sigma_all"].values)) + df["Mach_all"] = list(map(lambda v: v/cs0, df["sigma_all"].values)) + df["time"] = list(map(lambda x : x.T[0], + [pl.study.get_value(f"/series/time", unit=U.Myr)[run] + for run in pl.runs])) + + df["sigma"] = list(map(lambda l: np.mean(l), df["sigma_all"].values)) + df["Mach"] = df["sigma"] / cs0 + df["turnover"] = (df["L"] * U.pc.express(U.km) / (2 * df["sigma"]))* U.s.express(U.Myr) + + return df + + +def rho_pdf(pp): + pp.load_cells() + rho = pp.cells["rho"] * pp.info["unit_density"].express(U.H_cc) + rho_0 = np.mean(rho) + print(rho_0) + s = np.log(rho/rho_0) + values, edges = np.histogram(s, bins=300, range=(-15, 11), + density=True) + pp.unload_cells() + centers = 0.5 * (edges[1:] + edges[:-1]) + return (np.stack([values, centers]), {"logbins": True}) + +rule_pdf=Rule(rho_pdf, "Density PDF", name="rho_pdf", group="/hist") +def apply_rule_pdf(pp): + return pp.process(rule_pdf, pp, overwrite=True) diff --git a/utils/__init__.py b/utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/line_annotation.py b/utils/line_annotation.py similarity index 100% rename from line_annotation.py rename to utils/line_annotation.py diff --git a/mypool.py b/utils/mypool.py similarity index 100% rename from mypool.py rename to utils/mypool.py diff --git a/params.py b/utils/params.py similarity index 94% rename from params.py rename to utils/params.py index 99749db..3f7ca68 100644 --- a/params.py +++ b/utils/params.py @@ -33,7 +33,7 @@ def params_from_file(filename): def default_params(): - return params_from_file(_dir_path + "/params.yml") + return params_from_file(_dir_path + "/../params.yml") def load_params(filename): diff --git a/units.py b/utils/units.py similarity index 100% rename from units.py rename to utils/units.py