From 1a07edb5b70bfe769c1ecd6cca4dd043ca103364 Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Thu, 14 May 2020 15:55:01 +0200 Subject: [PATCH] Add preliminary power spectrum integration --- plotter.py | 21 ++++++++++ postprocessor.py | 7 ++++ pspec_new.py | 83 +++++++++++++++++++++++++++++++++---- pspec_read.py | 104 ++++++++++++++++++++++++++++++++--------------- 4 files changed, 175 insertions(+), 40 deletions(-) diff --git a/plotter.py b/plotter.py index 5130f82..3dd3849 100644 --- a/plotter.py +++ b/plotter.py @@ -17,6 +17,7 @@ if os.environ.get("DISPLAY", "") == "": mpl.use("Agg") import pylab as P from comparator import * +import pspec_read P.rcParams["image.cmap"] = "plasma" P.rcParams["savefig.dpi"] = 400 @@ -634,6 +635,8 @@ class Plotter(Aggregator, BaseProcessor): colors=None, nml_color=None, legend=None, + subname_x=None, + subname_y=None, **kwargs ): @@ -643,6 +646,13 @@ class Plotter(Aggregator, BaseProcessor): node_x = self.save.get_node(name_x) node_y = self.save.get_node(name_y) + if subname_x: + hdf5_x = tables.open_file(node_x.read()) + node_x = hdf5_x.get_node(subname_x) + if subname_y: + hdf5_y = tables.open_file(node_y.read()) + node_y = hdf5_y.get_node(subname_y) + xlabel, xunit_old, xunit = self._ax_label_unit( node_x, xlabel, xunit, xunit_coeff ) @@ -772,6 +782,16 @@ class Plotter(Aggregator, BaseProcessor): color=base_line.get_color(), label=fitlabel, ) + if subname_x: + hdf5_x.close() + if subname_y: + hdf5_y.close() + + def _pspec(self, name, **kwargs): + del kwargs["run"] + file_pspec = self.save.get_node("/hdf5/pspec").read() + num = self.save.root._v_attrs.num + getattr(pspec_read, "pspec_" + name)(file_pspec, ".", num, **kwargs) def _overlay_fit(self, x, y, yerr=None, kind="linear", label=None, **kwargs): if kind == "linear": @@ -943,6 +963,7 @@ class Plotter(Aggregator, BaseProcessor): partial(self._plot, "/profile/axis", "/profile/rho_prof"), dependencies=["axis", "rho_prof"], ), + "pspec": PlotRule(self, self._pspec, dependencies={"pspec": None}), } averageables = ["coldens", "rho", "T", "Q"] diff --git a/postprocessor.py b/postprocessor.py index ae8ef19..e40f5c3 100644 --- a/postprocessor.py +++ b/postprocessor.py @@ -1,5 +1,6 @@ # coding: utf-8 import pandas as pd +import pspec_new from baseprocessor import * mass_func = lambda dset: dset["rho"] * dset["dx"] ** 3 # Mass function @@ -584,6 +585,11 @@ class PostProcessor(HDF5Container): return sinks_dict + def _pspec(self): + outfile = self.path_out + "/pspec.h5" + pspec_new.pspec(repo=self.path, iouts=[self.num], outfile=outfile) + return outfile + def def_rules(self): self.rules = { @@ -674,6 +680,7 @@ class PostProcessor(HDF5Container): "Teff": cst.K, }, ), + "pspec": Rule(self, self._pspec, "Power spectrum", "/hdf5"), # Helpers "radial_bins": Rule(self, self._radial_bins, "Radial bins", "/radial"), "rr": Rule(self, self._rr, "Coordinate map", "/maps"), diff --git a/pspec_new.py b/pspec_new.py index 4ec8536..f48eaef 100644 --- a/pspec_new.py +++ b/pspec_new.py @@ -15,7 +15,6 @@ from pymses.analysis import cube3d import pymses.utils.misc import tables as T -import bunch from i_utils import args from i_utils import tools @@ -464,6 +463,45 @@ 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( + "repo", help="RAMSES output repository", type=str, default=".", nargs="?" +) +parser.add_argument( + "iouts", help="output numbers", type=args.selection, default=":", nargs="?" +) +parser.add_argument( + "outfile", + help="output file format (see below for fields)", + type=str, + default=".", + nargs="?", +) +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 ) @@ -504,7 +542,7 @@ def main(arg): parser.error("Invalid slicing axis: %r" % arg.sliceaxis) # Output numbers - iouts = tools.select_outputs(arg.repo, arg.iouts) + iouts = arg.iouts # tools.select_outputs(arg.repo, arg.iouts) for iout in iouts: print("Output %d" % iout) @@ -547,7 +585,7 @@ def main(arg): distance=arg.size / 2.0, far_cut_depth=arg.size / 2.0, up_vector="y", - map_max_size=256, + map_max_size=2 ** clvl, ) cubes = {} @@ -555,7 +593,7 @@ def main(arg): operator = ScalarOperator(cube_vars[i], cube_units[i]) extractor = cube3d.CubeExtractor(amr, operator) cubes[v] = extractor.process( - cam, cube_size=arg.size, resolution=256 + cam, cube_size=arg.size, resolution=2 ** clvl ).data else: h5f = T.open_file("cube.hdf5", "r") @@ -577,7 +615,6 @@ def main(arg): % (clvl, read_lvl) ) clvl = read_lvl - # Degrade cubes ---------------------------------------------------------------- if clvl < read_lvl: print("Degrade cubes") @@ -655,10 +692,23 @@ def main(arg): pcubes["krs"] = pcube(fcubes["krsx"], fcubes["krsy"], fcubes["krsz"]) pcubes["B"] = pcube(fcubes["Bx"], fcubes["By"], fcubes["Bz"]) pcubes["cos_vB"] = pcube(fcubes["cos_vB"]) + pcubes["vz"] = pcube(fcubes["vz"]) print("Compute 3D power spectra") # 3D power spectra ------------------------------------------------------------- - for v in ["rho", "logrho", "v", "kr", "vc", "vs", "krc", "krs", "B", "cos_vB"]: + for v in [ + "rho", + "logrho", + "v", + "kr", + "vc", + "vs", + "krc", + "krs", + "B", + "cos_vB", + "vz", + ]: pspec, kbins, pspec2, fbins = pspectrum( pcubes[v], cubes_k["k"], kbins, knorm, arg.fbins ) @@ -751,6 +801,7 @@ def main(arg): "By", "Bz", "cos_vB", + "vz", ]: fcubes2[v] = ifft(fcubes[v], axis=saxis) @@ -770,12 +821,25 @@ def main(arg): pcubes2["krs"] = pcube(fcubes2["krsx"], fcubes2["krsy"], fcubes2["krsz"]) pcubes2["B"] = pcube(fcubes2["Bx"], fcubes2["By"], fcubes2["Bz"]) pcubes2["cos_vB"] = pcube(fcubes2["cos_vB"]) + pcubes2["vz"] = pcube(fcubes2["vz"]) print("Compute 2D power spectra") # 2D power spectra ------------------------------------------------------------- ns = 2 ** clvl f = "_%%(i)0%dd" % (np.floor(np.log10(ns)) + 1) - for v in ["rho", "logrho", "v", "kr", "vc", "vs", "krc", "krs", "B", "cos_vB"]: + for v in [ + "rho", + "logrho", + "v", + "kr", + "vc", + "vs", + "krc", + "krs", + "B", + "cos_vB", + "vz", + ]: for i in xrange(ns): pspec, kbins, pspec2, fbins = pspectrum( pcubes2[v][:, :, i], cubes_k["k"][:, :, i], kbins, knorm, arg.fbins @@ -857,4 +921,7 @@ if __name__ == "__main__": def pspec(**kwargs): - main(bunch.bunchify(kwargs)) + arg = parser.parse_args("") + for kwarg in kwargs: + setattr(arg, kwarg, kwargs[kwarg]) + main(arg) diff --git a/pspec_read.py b/pspec_read.py index 7c00234..6268b9a 100644 --- a/pspec_read.py +++ b/pspec_read.py @@ -8,7 +8,7 @@ import matplotlib.pyplot as P ################################################################ -def pspec_rho(file, path_out, num, color="r", save=False): +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") @@ -22,13 +22,13 @@ def pspec_rho(file, path_out, num, color="r", save=False): 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) + 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): +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") @@ -42,13 +42,13 @@ def pspec_B(file, path_out, num, color="r", save=False): 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) + 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): +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") @@ -62,13 +62,33 @@ def pspec_v(file, path_out, num, color="r", save=False): 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) + 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_cos(file, path_out, num, color="r", save=False): +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") @@ -82,13 +102,13 @@ def pspec_cos(file, path_out, num, color="r", save=False): 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) + 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): +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") @@ -102,13 +122,13 @@ def pspec_logrho(file, path_out, num, color="r", save=False): 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) + 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): +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") @@ -122,13 +142,13 @@ def pspec_kr(file, path_out, num, color="r", save=False): 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) + 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): +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") @@ -142,13 +162,13 @@ def pspec_vc(file, path_out, num, color="r", save=False): 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) + 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): +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") @@ -162,7 +182,7 @@ def pspec_vs(file, path_out, num, color="r", save=False): 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) + 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") @@ -173,7 +193,7 @@ def pspec_vs(file, path_out, num, color="r", save=False): ################################################################### -def pspec_rho_2D(file, path_out, num, plan, color="r", save=False): +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") @@ -187,13 +207,13 @@ def pspec_rho_2D(file, path_out, num, plan, color="r", save=False): 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) + 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): +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") @@ -207,13 +227,13 @@ def pspec_B_2D(file, path_out, num, plan, color="r", save=False): 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) + 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): +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") @@ -227,13 +247,33 @@ def pspec_v_2D(file, path_out, num, plan, color="r", save=False): 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) + 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_cos_2D(file, path_out, num, plan, color="r", save=False): +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") @@ -247,13 +287,13 @@ def pspec_cos_2D(file, path_out, num, plan, color="r", save=False): 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) + 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): +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") @@ -267,13 +307,13 @@ def pspec_logrho_2D(file, path_out, num, plan, color="r", save=False): 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) + 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): +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") @@ -287,13 +327,13 @@ def pspec_kr_2D(file, path_out, num, plan, color="r", save=False): 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) + 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): +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") @@ -307,13 +347,13 @@ def pspec_vc_2D(file, path_out, num, plan, color="r", save=False): 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) + 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): +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") @@ -327,7 +367,7 @@ def pspec_vs_2D(file, path_out, num, plan, color="r", save=False): 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) + 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")