From 6b3f127878ca40d83a970901272737a25a7e921e Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Fri, 24 May 2019 16:29:27 +0200 Subject: [PATCH] Zoom --- disk_postprocess.py | 90 +++++++++++++++++++++++++++++++++------------ pipeline_disk.py | 52 +++++++++++++++++++++++++- 2 files changed, 117 insertions(+), 25 deletions(-) diff --git a/disk_postprocess.py b/disk_postprocess.py index 27f2085..bdc5bf7 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -35,6 +35,7 @@ P.rcParams["savefig.dpi"] = 400 def make_image_disk( path, num, + rad=0.5, path_out=None, order="<", save_data=True, @@ -67,7 +68,6 @@ def make_image_disk( ro = pymses.RamsesOutput(path, num, order=order) amr = ro.amr_source(["rho", "vel", "P"]) - rad = 0.5 center = [0.5, 0.5, 0.5] return make_image_aux( @@ -258,6 +258,7 @@ def make_image_aux( if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(name_im) P.close() @@ -323,6 +324,7 @@ def make_image_aux( if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(name_im) P.close() @@ -352,6 +354,7 @@ def make_image_aux( if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(name_im) P.close() @@ -398,7 +401,13 @@ def make_image_aux( # maps_disk['Q_' + ax_los] = map_Q im = P.imshow( - map_Q, extent=im_extent, origin="lower", norm=mpl.colors.LogNorm() + map_Q, + extent=im_extent, + origin="lower", + cmap="RdBu", + norm=mpl.colors.LogNorm(), + vmin=0.01, + vmax=100.0, ) P.locator_params(axis="x", nbins=ntick) @@ -416,6 +425,7 @@ def make_image_aux( if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(name_im) P.close() @@ -449,6 +459,7 @@ def make_image_aux( if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(name_im) P.close() @@ -696,7 +707,7 @@ def disk_prop( f.close() -def plot_disk_prop(path, num, force=False, tag="", interactive=False): +def plot_disk_prop(path, num, force=False, tag="", interactive=False, put_title=False): """ Plot properties of a disk @@ -740,6 +751,7 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False): if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + out_ext) P.close() @@ -749,10 +761,12 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False): P.plot(rad, prop_disk["temp"], color="k", linewidth=2) P.ylabel(r"$T \, (K)$") P.xlabel("disk radius") - P.title(title) + if put_title: + P.title(title) if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(path + "/T_disk_r_" + str(num).zfill(5) + out_ext) P.close() @@ -771,6 +785,7 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False): if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(path + "/V_disk_r_" + str(num).zfill(5) + out_ext) P.close() @@ -780,10 +795,12 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False): P.plot(rad, prop_disk["coldens"], color="k", linewidth=2) P.ylabel(r"$N\, (cm^{-2})$") P.xlabel("disk radius ") - P.title(title) + if put_title: + P.title(title) if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + out_ext) P.close() @@ -813,11 +830,12 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False): linewidth=2, label=r"$\alpha_{tot}$", ) - P.title( - title - + r", $\bar{\alpha}_{Reynolds} = %.1e, \bar{\alpha}_{grav} = %.1e$" - % (alpha_rey_mean, alpha_grav_mean) - ) + if put_title: + P.title( + title + + r", $\bar{\alpha}_{Reynolds} = %.1e, \bar{\alpha}_{grav} = %.1e$" + % (alpha_rey_mean, alpha_grav_mean) + ) P.legend() P.ylabel(r"$\alpha$") @@ -826,23 +844,26 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False): if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + out_ext) P.close() # Q - P.ylim([0, 10.0]) - P.xlim([0, 0.5]) - P.yticks(np.arange(0.0, 11, 1.0)) + P.ylim([0, 3.0]) + P.xlim([0, 0.25]) + P.yticks(np.arange(0.0, 3, 1.0)) P.grid() P.plot(rad, abs(prop_disk["Q_kepl"]), color="b", linewidth=2) - P.plot(rad, abs(prop_disk["Q_mean"]) * np.ones(len(rad)), "b:", linewidth=1) + # P.plot(rad, abs(prop_disk['Q_mean']) * np.ones(len(rad)), 'b:', linewidth=1) P.ylabel(r"$Q$") P.xlabel("disk radius ") - P.title(title + ", mass of disk = {} (code)".format(mass)) + if put_title: + P.title(title + ", mass of disk = {} (code)".format(mass)) if interactive: pass else: + P.tight_layout(pad=1) P.savefig(path + "/Q_r_" + str(num).zfill(5) + out_ext) P.close() @@ -851,11 +872,13 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False): P.plot(rad, abs(prop_disk["height"] / rad), color="b", linewidth=2) P.ylabel(r"H ratio") P.xlabel("disk radius ") - P.title(title + ", mass of box = {} (code)".format(prop_disk["mass_box"])) + if put_title: + P.title(title + ", mass of box = {} (code)".format(prop_disk["mass_box"])) if interactive: pass else: + P.tight_layout(pad=1) P.savefig(path + "/H_r_" + str(num).zfill(5) + out_ext) P.close() @@ -958,7 +981,7 @@ def disk_pdf( P.yscale("log") P.ylim([0.5 / nb_cells, 1.0]) P.xlabel(r"$\log(N / \bar{N})$") - P.ylabel(r"\mathcal{P}_N") + P.ylabel(r"$\mathcal{P}_N$") if put_title: P.title(title) values, edges, _ = P.hist( @@ -972,14 +995,17 @@ def disk_pdf( # Compute the slope of the right part of the histogramm mask_fit = (centers > 0.25) & (centers < 1.25) & (values > 0) if np.sum(mask_fit > 0): - (a, b, rho, _, _) = linregress(centers[mask_fit], np.log10(values[mask_fit])) - P.plot(centers, 10 ** (a * centers + b), "--") + (a, b, rho, _, stderr) = linregress( + centers[mask_fit], np.log10(values[mask_fit]) + ) + P.plot(centers, 10 ** (a * centers + b), "--", linewidth=2) print("a=%e, b=%e, rho=%e" % (a, b, rho)) fit = { "beta": int(tag.split("_")[1][4:]), "slope": a, "origin": b, "correlation": rho, + "stderr": stderr, } f = open(name_prop, "w") prop_disk["fit"] = fit @@ -989,6 +1015,7 @@ def disk_pdf( if interactive: pass else: + P.tight_layout(pad=1) P.savefig(path + "/dcol_hist_" + tag + "_" + str(num).zfill(5) + out_ext) P.close() @@ -1010,6 +1037,7 @@ def disk_pdf( if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(name_im) P.close() @@ -1049,11 +1077,14 @@ def compare( if pdf: beta = np.zeros(len(runs)) slope = np.zeros(len(runs)) + slope_std = np.zeros(len(runs)) for i, run in enumerate(runs): path_run = path + "/" + run nb_outputs = 0 + slope_run = [] + for num in nums: try: # Load property file @@ -1076,7 +1107,7 @@ def compare( if pdf: fit = prop_disk["fit"] beta[i] = fit["beta"] - slope[i] = slope[i] + fit["slope"] + slope_run.append(fit["slope"]) nb_outputs = nb_outputs + 1 print(run, num, nb_outputs) @@ -1089,7 +1120,8 @@ def compare( alpha_grav[i] = alpha_grav[i] / nb_outputs Q[i] = Q[i] / nb_outputs if pdf: - slope[i] = slope[i] / nb_outputs + slope[i] = np.mean(slope_run) + slope_std[i] = np.std(slope_run) else: for array in [alpha_rey, alpha_grav, Q]: array[i] = np.nan @@ -1119,6 +1151,7 @@ def compare( if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(path_out + "/alphaQ_" + nums_name + out_ext) P.close() @@ -1138,21 +1171,27 @@ def compare( if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(path_out + "/alphaQ0_" + nums_name + out_ext) P.close() if pdf: # slope of the pdf = f(beta) P.grid() - P.plot(beta, slope, "o-.") + P.errorbar(beta, slope, yerr=slope_std, fmt="o") P.legend() - P.ylabel(r"$d\log\% / d\logN$") + P.ylabel(r"$d\log\mathcal{P}_N / d\logN$") P.xlabel(r"$\beta$") + (a, b, rho, _, stderr) = linregress(beta, slope) + P.plot(beta, a * beta + b, "--", linewidth=1.5) + print("a=%e, b=%e, rho^2=%e" % (a, b, rho ** 2)) + if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(path_out + "/dcolslopebeta_" + nums_name + out_ext) P.close() @@ -1175,6 +1214,7 @@ def compare( if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(path_out + "/alphabeta_" + nums_name + out_ext) P.close() @@ -1248,6 +1288,7 @@ def evolution(path, nums, force=False, interactive=False, pdf=False): if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(path + "/alpha_time" + out_ext) P.close() @@ -1263,6 +1304,7 @@ def evolution(path, nums, force=False, interactive=False, pdf=False): if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(path + "/Q_time" + out_ext) P.close() @@ -1278,6 +1320,7 @@ def evolution(path, nums, force=False, interactive=False, pdf=False): if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(path + "/mass_time" + out_ext) P.close() @@ -1292,5 +1335,6 @@ def evolution(path, nums, force=False, interactive=False, pdf=False): if interactive: P.figure() else: + P.tight_layout(pad=1) P.savefig(path + "/dcolslope_time" + out_ext) P.close() diff --git a/pipeline_disk.py b/pipeline_disk.py index 704c638..95caa80 100644 --- a/pipeline_disk.py +++ b/pipeline_disk.py @@ -44,6 +44,40 @@ parser.add_argument( parser.add_argument("-i", "--interactive", help="Interactive mode", action="store_true") +parser.add_argument( + "-d", "--disk", help="compute specific disk radial analysis", action="store_true" +) +parser.add_argument( + "-pd", "--plot_disk", help="plot specific disk radial analysis", action="store_true" +) +parser.add_argument("-m", "--maps", help="do generic maps", action="store_true") +parser.add_argument( + "-c", "--compare", help="compare different runs", action="store_true" +) +parser.add_argument( + "-ev", + "--evolution", + help="plot evolution of quantities over time", + action="store_true", +) +parser.add_argument( + "--pdf", help="plot pdf of fluctuations of column density", action="store_true" +) + +parser.add_argument( + "--fft", help="use quick and dirty fft rendering", action="store_true" +) +parser.add_argument("--level", help="plot levels", action="store_true") +parser.add_argument("--cpu", help="plot cpu", action="store_true") +parser.add_argument("--zoom", help="zoom", type=float, default=2.0) +parser.add_argument( + "-ms", + "--mapsize", + help="size of the maps created in he map mode (in pixel)", + type=int, + default=1024, +) + parser.add_argument( "-d", "--disk", help="do specific disk radial analysis", action="store_true" ) @@ -108,10 +142,11 @@ parser.add_argument( parser.add_argument( "--format", help="Format of the output", - choices=["png", "jpeg", "pdf", "ps"], + choices=["png", "jpeg", "pdf", "svg", "ps"], default="jpeg", ) parser.add_argument("--dpi", help="Resolution of the output", type=int, default=400) +parser.add_argument("--beamer", help="Beamer mode", action="store_true") args = parser.parse_args() @@ -123,8 +158,17 @@ first = args.begin last = args.end step = args.step +rad = 0.5 / args.zoom + # extension for out files dp.out_ext = "." + args.format +if format == "pdf": + dp.P.style.use("pdf") + +if args.beamer: + dp.P.rcParams["font.family"] = "sans-serif" + dp.P.rcParams["figure.figsize"] = (5, 3.5) + # Plot properties dp.P.rcParams["image.cmap"] = args.colormap dp.P.rcParams["savefig.dpi"] = args.dpi @@ -159,6 +203,7 @@ for run in runs: maps_disk = dp.make_image_disk( path_in, i, + rad=rad, path_out=path_out, tag=run, map_size=args.mapsize, @@ -168,6 +213,7 @@ for run in runs: pos_star=np.array([args.x, args.y, args.z]), fft=args.fft, interactive=args.interactive, + put_title=(not args.beamer), ) print("[{}, {}] maps computed".format(run, i)) if args.disk: @@ -183,7 +229,7 @@ for run in runs: pos_star=np.array([args.x, args.y, args.z]), ) print("[{}, {}] disk_props computed".format(run, i)) - if args.disk: + if args.plot_disk: print("[{}, {}] plotting disk props".format(run, i)) dp.plot_disk_prop( path_out, @@ -191,6 +237,7 @@ for run in runs: tag=run, force=args.force_redo, interactive=args.interactive, + put_title=(not args.beamer), ) print("[{}, {}] disk props plotted".format(run, i)) if args.pdf: @@ -203,6 +250,7 @@ for run in runs: force=args.force_redo, tag=run, interactive=args.interactive, + put_title=(not args.beamer), ) print("[{}, {}] pdf computed".format(run, i)) # If we are here, success !