From 5eb5ae36c4dadf4ddd5f2c5bca49f6d491e1be4d Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Wed, 15 May 2019 11:41:29 +0200 Subject: [PATCH] linear interpolation + regression --- disk_postprocess.py | 131 ++++++++++++++++++++++++++++++++------------ pipeline_disk.py | 10 +++- 2 files changed, 106 insertions(+), 35 deletions(-) diff --git a/disk_postprocess.py b/disk_postprocess.py index 191362f..9387626 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -19,7 +19,7 @@ except: print("cPickle not found, using pickle") import pickle import tables - +from scipy.stats import linregress from pymses.sources.ramses import output from pymses.analysis import Camera, raytracing, slicing, splatting from pymses.filters import CellsToPoints @@ -872,8 +872,9 @@ def disk_pdf( interactive=False, nb_bin_hist=50, tag="", - rad_min=0.08, + rad_min=0.05, rad_max=0.4, + put_title=True, ): # Load property file @@ -900,7 +901,7 @@ def disk_pdf( print("maps file loaded") time = prop_disk["time"] - title = "t=" + str(time)[0:5] + " (code)" + title = tag.split("_")[1] + " t=" + str(time)[0:5] + " (code)" coldens_map = maps_disk["coldens_z"] im_extent = maps_disk["im_extent"] @@ -945,39 +946,58 @@ def disk_pdf( # Histogramm : density fluctuation distribution function dcoldens = np.log10(coldens_map.flatten() / coldens_mean) + mask_hist = (rr_flat > rad_min) & (rr_flat < rad_max) + nb_cells = np.sum(mask_hist) P.grid() P.yscale("log") + P.ylim([0.5 / nb_cells, 1.0]) P.xlabel(r"$\log(N / \bar{N})$") P.ylabel("Probability density") - P.title(title) - P.hist( - dcoldens[(rr_flat > rad_min) & (rr_flat < rad_max)], nb_bin_hist, range=(-1, 3) + if put_title: + P.title(title) + values, edges, _ = P.hist( + dcoldens[mask_hist], + nb_bin_hist, + range=(-1, 3), + weights=np.ones(nb_cells) / nb_cells, ) + centers = 0.5 * (edges[1:] + edges[:-1]) + + # Compute the slope of the right part of the histogramm + mask_fit = (centers > 0.25) & (centers < 1.25) & (values > 0) + (a, b, rho, _, _) = linregress(centers[mask_fit], np.log10(values[mask_fit])) + print("a=%e, b=%e, rho=%e" % (a, b, rho)) + fit = { + "beta": int(tag.split("_")[1][4:]), + "slope": a, + "origin": b, + "correlation": rho, + } + P.plot(centers, 10 ** (a * centers + b), "--") + f = open(name_prop, "w") + prop_disk["fit"] = fit + pickle.dump(prop_disk, f) + f.close() if interactive: pass else: - P.savefig(path + "/dcol_hist_" + str(num).zfill(5) + out_ext) + P.savefig(path + "/dcol_hist_" + tag + "_" + str(num).zfill(5) + out_ext) P.close() # Fluctuation map dcoldens_map = np.reshape(dcoldens, coldens_map.shape) - # cont = P.contour(coldens_mean_map, - # extent=im_extent, - # origin='lower', - # colors='k', - # linewidths=0.1) - im = P.imshow(dcoldens_map, extent=im_extent, origin="lower", cmap="viridis") - P.title(title) + if put_title: + P.title(title) P.xlabel(r"$x$") P.ylabel(r"$y$") cbar = P.colorbar(im) cbar.set_label(r"$log(N/\bar{N})$ (code)") - name = path + "/dcoldensmap_" + format(num, "05") + name = path + "/dcoldensmap_" + tag + "_" + format(num, "05") name_im = name + out_ext if interactive: @@ -987,7 +1007,7 @@ def disk_pdf( P.close() -def compare(path, runs, num, force=False, interactive=False): +def compare(path, runs, num, force=False, interactive=False, Q_in_name=True, pdf=False): """ Compare properties of a disk in several simulations @@ -1003,7 +1023,11 @@ def compare(path, runs, num, force=False, interactive=False): alpha_rey = np.zeros(len(runs)) alpha_grav = np.zeros(len(runs)) Q = np.zeros(len(runs)) - Q0 = np.zeros(len(runs)) + if Q_in_name: + Q0 = np.zeros(len(runs)) + if pdf: + beta = np.zeros(len(runs)) + slope = np.zeros(len(runs)) for i, run in enumerate(runs): path_run = path + "/" + run @@ -1021,12 +1045,17 @@ def compare(path, runs, num, force=False, interactive=False): alpha_rey[i] = prop_disk["alpha_rey_mean"] alpha_grav[i] = prop_disk["alpha_grav_mean"] Q[i] = prop_disk["Q_min"] - Q0[i] = float(run.split("_")[2][1:]) + if Q_in_name: + Q0[i] = float(run.split("_")[2][1:]) + if pdf: + fit = prop_disk["fit"] + beta[i] = fit["beta"] + slope[i] = fit["slope"] # Check if the output file exists, and exit if it is the case name_save = path + "/alphaQ_" + str(num).zfill(5) + out_ext - if not force and len(glob.glob(name_save)) != 0: - return + # if (not force and len(glob.glob(name_save)) !=0): + # return title = "t=" + str(time[0]) + " (code)" @@ -1048,23 +1077,57 @@ def compare(path, runs, num, force=False, interactive=False): P.savefig(path + "/alphaQ_" + str(num).zfill(5) + out_ext) P.close() - # alpha = f(Q0) - P.yscale("log") - P.ylim([1e-7, 1.0]) - P.grid() + if Q_in_name: + # alpha = f(Q0) + P.yscale("log") + P.ylim([1e-7, 1.0]) + P.grid() - P.plot(Q0, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$") - P.plot(Q0, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$") + P.plot(Q0, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$") + P.plot(Q0, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$") - P.legend() - P.ylabel(r"$\bar{\alpha}$") - P.xlabel(r"$Q_{0}$") + P.legend() + P.ylabel(r"$\bar{\alpha}$") + P.xlabel(r"$Q_{0}$") - if interactive: - P.figure() - else: - P.savefig(path + "/alphaQ0_" + str(num).zfill(5) + out_ext) - P.close() + if interactive: + P.figure() + else: + P.savefig(path + "/alphaQ0_" + str(num).zfill(5) + out_ext) + P.close() + + if pdf: + # slope of the pdf = f(beta) + P.grid() + P.plot(beta, slope, "o-.") + + P.legend() + P.ylabel(r"$d\log\% / d\logN$") + P.xlabel(r"$\beta$") + + if interactive: + P.figure() + else: + P.savefig(path + "/dcolslopebeta_" + str(num).zfill(5) + out_ext) + P.close() + + # alpha = f(beta) + P.yscale("log") + P.ylim([1e-7, 1.0]) + P.grid() + + P.plot(beta, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$") + P.plot(beta, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$") + + P.legend() + P.ylabel(r"$\bar{\alpha}$") + P.xlabel(r"$\beta_{cool}$") + + if interactive: + P.figure() + else: + P.savefig(path + "/alphabeta_" + str(num).zfill(5) + out_ext) + P.close() def evolution(path, nums, force=False, interactive=False): diff --git a/pipeline_disk.py b/pipeline_disk.py index 68b20b6..3f241b2 100644 --- a/pipeline_disk.py +++ b/pipeline_disk.py @@ -224,4 +224,12 @@ if args.compare: path = storage_out + path_suffix for i in range(first, last + 1, step): - dp.compare(path, runs, i, force=args.force_redo, interactive=args.interactive) + dp.compare( + path, + runs, + i, + force=args.force_redo, + interactive=args.interactive, + Q_in_name=(not args.pdf), + pdf=args.pdf, + )