Zoom

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")
+    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 ")
+    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,6 +830,7 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
         linewidth=2,
         label=r"$\alpha_{tot}$",
     )
+    if put_title:
         P.title(
             title
             + r", $\bar{\alpha}_{Reynolds} = %.1e,  \bar{\alpha}_{grav} = %.1e$"
@@ -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 ")
+    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 ")
+    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()

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 !