Added comparaison and evolution mode

disk_postprocess.py

@@ -549,6 +549,9 @@ def disk_prop(
     alpha_grav_rad = np.zeros(nb_bin - 1)
     Q_kepl_rad = np.zeros(nb_bin - 1)
     height_rad = np.zeros(nb_bin - 1)
+    vol_rad = np.zeros(nb_bin - 1)  # Volume of a bin
+    surf_rad = np.zeros(nb_bin - 1)  # Surface of a bin
+    mass_rad = np.zeros(nb_bin - 1)  # Mass of a bin
 
     # Density fluctuations
     hist_drho = np.zeros(nb_bin_hist)
@@ -559,34 +562,33 @@ def disk_prop(
 
         # print("Bin #{} : {} cells between {} and {}".format(i, np.sum(mask_bin), rad[i], rad[i + 1]))
 
-        press_rad[i] = np.sum(press_disk[mask_bin] * dvol_disk[mask_bin]) / np.sum(
+        vol_rad[i] = np.sum(dvol_disk[mask_bin])
-            dvol_disk[mask_bin]
+        mass_rad[i] = np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin])
-        )
+        press_rad[i] = np.sum(press_disk[mask_bin] * dvol_disk[mask_bin]) / vol_rad[i]
-        rho_rad[i] = np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) / np.sum(
+        rho_rad[i] = np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) / vol_rad[i]
-            dvol_disk[mask_bin]
+        temp_rad[i] = np.sum(press_disk[mask_bin] * dvol_disk[mask_bin]) / mass_rad[i]
-        )
-        temp_rad[i] = np.sum(press_disk[mask_bin] * dvol_disk[mask_bin]) / np.sum(
-            rho_disk[mask_bin] * dvol_disk[mask_bin]
-        )
 
         # Surface of a bin :  S = dr * 2 * pi * r with
         #  dr = rad[i + 1] - rad[i] and r = (rad[i + 1] + rad[i]) / 2.
-        coldens_rad[i] = np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) / (
+        surf_rad[i] = (rad[i + 1] - rad[i]) * (rad[i + 1] + rad[i]) * np.pi
-            (rad[i + 1] - rad[i]) * (rad[i + 1] + rad[i]) * np.pi
+        coldens_rad[i] = mass_rad[i] / surf_rad[i]
+
+        v_az_rad[i] = (
+            np.sum(v_az_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin])
+            / mass_rad[i]
         )
 
-        v_az_rad[i] = np.sum(
+        v_rad_rad[i] = (
-            v_az_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin]
+            np.sum(v_rad_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin])
-        ) / np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin])
+            / mass_rad[i]
+        )
 
-        v_rad_rad[i] = np.sum(
+        height_rad[i] = (
-            v_rad_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin]
+            np.sum(height_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin])
-        ) / np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin])
+            / mass_rad[i]
-        height_rad[i] = np.sum(
+        )
-            height_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin]
-        ) / np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin])
 
-        alpha_rey_rad[i] = (
+        alpha_rey_rad[i] = (2.0 / 3) * (
             (
                 np.sum(
                     v_az_disk[mask_bin]
@@ -601,12 +603,6 @@ def disk_prop(
             / abs(v_az_rad[i])
         )
 
-        # alpha_rey_rad_bis[i] = (2./3) * (np.sum((v_az_disk[mask_bin] - v_az_rad[i])
-        #                                         * (v_rad_disk[mask_bin] - v_rad_rad[i])
-        #                                         * rho_disk[mask_bin] * dvol_disk[mask_bin])
-        #                                  / np.sum(dvol_disk[mask_bin] * press_disk[mask_bin])
-        #                                  * v_az_rad[i] / abs(v_az_rad[i]))
-
         alpha_grav_rad[i] = (2.0 / 3) * (
             np.sum(
                 g_az_disk[mask_bin]
@@ -619,9 +615,10 @@ def disk_prop(
             * coldens_rad[i]
         )
 
-        v_kepl_rad[i] = np.sum(
+        v_kepl_rad[i] = (
-            v_kepl[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin]
+            np.sum(v_kepl[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin])
-        ) / np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin])
+            / mass_rad[i]
+        )
 
         # Histogramm : density fluctuaction distribution function
         drho = np.log(rho_disk[mask_bin] / rho_rad[i])
@@ -630,19 +627,33 @@ def disk_prop(
         )
         hist_drho = hist_drho + hist
 
-    print(hist_drho, hist_edges)
+    # Derived quantities
     cs_rad = np.sqrt(temp_rad)
     Q_kepl_rad = cs_rad * v_az_rad / (np.pi * G * coldens_rad * rad[0 : nb_bin - 1])
 
+    # Means
+    mask_mean = (0.1 < rad[0 : nb_bin - 1]) & (rad[0 : nb_bin - 1] < 0.2)
+    print(rad[0 : nb_bin - 1][mask_mean])
+    mass_mean = np.sum(mass_rad[mask_mean])
+    alpha_rey_mean = np.sum(alpha_rey_rad[mask_mean] * mass_rad[mask_mean]) / mass_mean
+    alpha_grav_mean = (
+        np.sum(alpha_grav_rad[mask_mean] * mass_rad[mask_mean]) / mass_mean
+    )
+    Q_mean = np.sum(Q_kepl_rad[mask_mean] * mass_rad[mask_mean]) / mass_mean
+    Q_min = np.nanmin(Q_kepl_rad)
+    print("alphas, Q ", alpha_rey_mean, alpha_grav_mean, Q_mean)
+
+    # store the results
     prop_disk = {
         "time": time,
-        "tot_mass": total_mass_disk,
+        "mass_disk": total_mass_disk,
         "mass_box": total_mass,
         "rad": rad[0 : nb_bin - 1],
         "center": pos_star,
         "alpha_rey": alpha_rey_rad,
-        # 'alpha_rey_bis':alpha_rey_rad_bis,
+        "alpha_rey_mean": alpha_rey_mean,
         "alpha_grav": alpha_grav_rad,
+        "alpha_grav_mean": alpha_grav_mean,
         "v_rad": v_rad_rad,
         "v_az": v_az_rad,
         "v_kepl": v_kepl_rad,
@@ -654,10 +665,10 @@ def disk_prop(
         "temp": temp_rad,
         "cs": cs_rad,
         "Q_kepl": Q_kepl_rad,
+        "Q_mean": Q_mean,
+        "Q_min": Q_min,
         "height": height_rad,
     }
-
-    # store the results
     f = open(name_save, "w")
     pickle.dump(prop_disk, f)
     f.close()
@@ -677,7 +688,7 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
 
     # Check if the properties file exists
     if len(glob.glob(name_save)) == 0:
-        raise ("no pickle file for disk properties. Run single_z_disk_prop")
+        raise ("no pickle file for disk properties. Run disk_prop() first")
     f = open(name_save, "r")
     prop_disk = pickle.load(f)
     f.close()
@@ -688,8 +699,9 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
         return
 
     time = prop_disk["time"]
-    mass = prop_disk["tot_mass"]
+    mass = prop_disk["mass_disk"]
     title = "t=" + str(time)[0:5] + " (code)"
+    rad = prop_disk["rad"]
 
     if interactive:
         P.figure()
@@ -699,7 +711,7 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
     P.xscale("log")
     P.yscale("log")
     P.grid()
-    P.plot(prop_disk["rad"], prop_disk["rho"], color="k", linewidth=2)
+    P.plot(rad, prop_disk["rho"], color="k", linewidth=2)
     P.ylabel(r"$n \, (code)$")
     P.xlabel("disk radius")
     P.title(title)
@@ -712,7 +724,7 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
     P.xscale("log")
     P.yscale("log")
     P.grid()
-    P.plot(prop_disk["rad"], prop_disk["temp"], color="k", linewidth=2)
+    P.plot(rad, prop_disk["temp"], color="k", linewidth=2)
     P.ylabel(r"$T \, (K)$")
     P.xlabel("disk radius")
     P.title(title)
@@ -725,10 +737,10 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
     P.xscale("log")
     P.yscale("symlog", linthreshy=0.01)
 
-    P.plot((prop_disk["rad"]), ((prop_disk["v_rad"])), color="k", linewidth=2)
+    P.plot(rad, prop_disk["v_rad"], color="k", linewidth=2)
-    P.plot((prop_disk["rad"]), ((prop_disk["v_kepl"])), color="b", linewidth=2)
+    P.plot(rad, prop_disk["v_kepl"], color="b", linewidth=2)
-    P.plot((prop_disk["rad"]), (abs(prop_disk["v_az"])), color="r", linewidth=2)
+    P.plot(rad, abs(prop_disk["v_az"]), color="r", linewidth=2)
-    P.plot((prop_disk["rad"]), ((prop_disk["cs"])), color="c", linewidth=2)
+    P.plot(rad, prop_disk["cs"], color="c", linewidth=2)
     P.grid()
     P.legend((r"$v_r$", r"$v_{kepl}$", r"$v_\phi$", r"$c_s$"), loc="upper right")
 
@@ -743,7 +755,7 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
     P.xscale("log")
     P.yscale("log")
     P.grid()
-    P.plot(prop_disk["rad"], prop_disk["coldens"], color="k", linewidth=2)
+    P.plot(rad, prop_disk["coldens"], color="k", linewidth=2)
     P.ylabel(r"$N\, (cm^{-2})$")
     P.xlabel("disk radius ")
     P.title(title)
@@ -754,36 +766,40 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
         P.close()
 
     # Alpha
+    alpha_rey_mean, alpha_grav_mean = (
+        prop_disk["alpha_rey_mean"],
+        prop_disk["alpha_grav_mean"],
+    )
+
     # P.xscale('log')
     P.xlim([1e-2, 0.25])
     P.yscale("log")
     P.ylim([1e-7, 1.0])
     P.grid()
     P.plot(
-        prop_disk["rad"],
+        rad, abs(prop_disk["alpha_rey"]), "b", linewidth=2, label=r"$\alpha_{Reynolds}$"
-        abs(prop_disk["alpha_rey"]),
-        linewidth=2,
-        label=r"$\alpha_{Reynolds}$",
     )
-    # P.plot(prop_disk['rad'],abs(prop_disk['alpha_rey_bis']), '--', linewidth=1,label=r"$\alpha_R 2$")
+    P.plot(rad, abs(alpha_rey_mean * np.ones(len(rad))), "b:", linewidth=1)
     P.plot(
-        prop_disk["rad"],
+        rad, abs(prop_disk["alpha_grav"]), "r", linewidth=2, label=r"$\alpha_{grav}$"
-        abs(prop_disk["alpha_grav"]),
-        linewidth=2,
-        label=r"$\alpha_{grav}$",
     )
+    P.plot(rad, abs(alpha_grav_mean * np.ones(len(rad))), "r:", linewidth=1)
     P.plot(
-        prop_disk["rad"],
+        rad,
         abs(prop_disk["alpha_rey"]) + abs(prop_disk["alpha_grav"]),
-        "--",
+        "g--",
         linewidth=2,
         label=r"$\alpha_{tot}$",
     )
+    P.title(
+        title
+        + r", $\bar{\alpha}_{Reynolds} = %.1e,  \bar{\alpha}_{grav} = %.1e$"
+        % (alpha_rey_mean, alpha_grav_mean)
+    )
 
     P.legend()
     P.ylabel(r"$\alpha$")
     P.xlabel("disk radius ")
-    P.title(title)
 
     if interactive:
         P.figure()
@@ -796,7 +812,8 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
     P.xlim([0, 0.5])
     P.yticks(np.arange(0.0, 11, 1.0))
     P.grid()
-    P.plot(prop_disk["rad"], abs(prop_disk["Q_kepl"]), color="b", linewidth=2)
+    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.ylabel(r"$Q$")
     P.xlabel("disk radius ")
     P.title(title + ", mass of disk = {} (code)".format(mass))
@@ -809,12 +826,7 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
 
     # height ratio
     P.grid()
-    P.plot(
+    P.plot(rad, abs(prop_disk["height"] / rad), color="b", linewidth=2)
-        prop_disk["rad"],
-        abs(prop_disk["height"] / prop_disk["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"]))
@@ -841,3 +853,176 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
     else:
         P.savefig(path + "/drho_hist_" + str(num).zfill(5) + out_ext)
         P.close()
+
+
+def compare(path, runs, num, force=False, interactive=False):
+    """
+    Compare properties of a disk in several simulations
+
+    num         id of the ramses output
+    runs        list of runs to consider
+    path        path to the properties file
+    force       if set, redo plots even if already done
+    interactive interactive mode, to use in a %pylab ipython shell
+    """
+
+    # Initialize arrays
+    time = np.zeros(len(runs))
+    alpha_rey = np.zeros(len(runs))
+    alpha_grav = np.zeros(len(runs))
+    Q = np.zeros(len(runs))
+    Q0 = np.zeros(len(runs))
+
+    for i, run in enumerate(runs):
+        path_run = path + "/" + run
+        # Load property file
+        name_save = path_run + "/prop_disk_" + str(num).zfill(5) + ".save"
+
+        # Check if the properties file exists
+        if len(glob.glob(name_save)) == 0:
+            raise ("no pickle file for disk properties. Run disk_prop() first")
+        f = open(name_save, "r")
+        prop_disk = pickle.load(f)
+        f.close()
+
+        time[i] = prop_disk["time"]
+        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:])
+
+    # 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
+
+    title = "t=" + str(time[0]) + " (code)"
+
+    # alpha = f(Qmin)
+    P.yscale("log")
+    P.ylim([1e-7, 1.0])
+    P.grid()
+
+    P.plot(Q, abs(alpha_rey), "o--", label=r"$\bar{\alpha}_{Reynolds}$")
+    P.plot(Q, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$")
+
+    P.legend()
+    P.ylabel(r"$\bar{\alpha}$")
+    P.xlabel(r"$Q_{min}$")
+
+    if interactive:
+        P.figure()
+    else:
+        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()
+
+    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}$")
+
+    if interactive:
+        P.figure()
+    else:
+        P.savefig(path + "/alphaQ0_" + str(num).zfill(5) + out_ext)
+        P.close()
+
+
+def evolution(path, nums, force=False, interactive=False):
+    """
+    Plot properties over time
+
+    path        path to the properties file
+    nums        list of id of the ramses output
+    force       if set, redo plots even if already done
+    interactive interactive mode, to use in a %pylab ipython shell
+    """
+
+    # Initialize arrays
+    time = np.zeros(len(nums))
+    alpha_rey = np.zeros(len(nums))
+    alpha_grav = np.zeros(len(nums))
+    Qmin = np.zeros(len(nums))
+    Qmean = np.zeros(len(nums))
+    mass_disk = np.zeros(len(nums))
+    mass_box = np.zeros(len(nums))
+
+    for i, num in enumerate(nums):
+
+        # Load property file
+        name_prop = path + "/prop_disk_" + str(num).zfill(5) + ".save"
+
+        # Check if the properties file exists
+        if len(glob.glob(name_prop)) == 0:
+            raise ("no pickle file for disk properties. Run disk_prop() first")
+        f = open(name_prop, "r")
+        prop_disk = pickle.load(f)
+        f.close()
+
+        time[i] = prop_disk["time"]
+        alpha_rey[i] = prop_disk["alpha_rey_mean"]
+        alpha_grav[i] = prop_disk["alpha_grav_mean"]
+        Qmin[i] = prop_disk["Q_min"]
+        Qmean[i] = prop_disk["Q_mean"]
+        mass_disk[i] = prop_disk["mass_disk"]
+        mass_box[i] = prop_disk["mass_box"]
+
+    # Check if the output file exists, and exit if it is the case
+    name_save = path + "/alpha_time" + out_ext
+    if not force and len(glob.glob(name_save)) != 0:
+        return
+
+    # Alpha
+    P.yscale("log")
+    P.ylim([1e-7, 1.0])
+    P.grid()
+
+    P.plot(time, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$")
+    P.plot(time, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$")
+
+    P.legend()
+    P.ylabel(r"$\bar{\alpha}$")
+    P.xlabel(r"time (code)")
+
+    if interactive:
+        P.figure()
+    else:
+        P.savefig(path + "/alpha_time" + out_ext)
+        P.close()
+
+    # Q
+    P.grid()
+    P.plot(time, Qmin, "o-.", label=r"$Q_{min}$")
+    P.plot(time, Qmean, "*--", label=r"$\bar{Q}$")
+
+    P.legend()
+    P.ylabel(r"$Q$")
+    P.xlabel(r"time (code)")
+
+    if interactive:
+        P.figure()
+    else:
+        P.savefig(path + "/Q_time" + out_ext)
+        P.close()
+
+    # M
+    P.grid()
+    P.plot(time, mass_disk, "o-.", label=r"$M_{disk}$")
+    P.plot(time, mass_box, "*--", label=r"$M_{box}$")
+
+    P.legend()
+    P.ylabel(r"$M / M_{*}$")
+    P.xlabel(r"time (code)")
+
+    if interactive:
+        P.figure()
+    else:
+        P.savefig(path + "/mass_time" + out_ext)
+        P.close()

pipeline_disk.py

@@ -12,35 +12,22 @@ storage_in = "/home/nbrucy/simus/"
 storage_out = "/home/nbrucy/visus/"
 
 parser = argparse.ArgumentParser()
+
 parser.add_argument("runs", help="name of runs", nargs="*", default=["015_iso"])
 parser.add_argument("-b", "--begin", help="id of first output", type=int, default=1)
 parser.add_argument("-e", "--end", help="id of last output", type=int, default=100)
 parser.add_argument("-s", "--step", help="step between two output", type=int, default=1)
-parser.add_argument(
-    "-d", "--disk", help="do specific disk radial analysis", action="store_true"
-)
-parser.add_argument(
-    "--nb_bin", help="Number of bins for azimuthal averages", type=int, default=50
-)
-
-
-parser.add_argument("-m", "--maps", help="do generic maps", action="store_true")
-
 parser.add_argument("-p", "--project", help="specify project name", default="disk")
-
 parser.add_argument(
     "-fr",
     "--force_redo",
     help="redo plots even if the files already exist",
     action="store_true",
 )
-
 parser.add_argument(
     "-w", "--watch", help="wait and watch for missing outputs", action="store_true"
 )
-
 parser.add_argument("--skip", help="skip failed loadings", action="store_true")
-
 parser.add_argument(
     "-wt",
     "--waiting_time",
@@ -54,9 +41,33 @@ parser.add_argument(
     help="number of allowed failures when waiting",
     default=30,
 )
+parser.add_argument("-i", "--interactive", help="Interactive mode", action="store_true")
+
+
+parser.add_argument(
+    "-d", "--disk", help="do 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(
+    "--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(
+    "--nb_bin", help="Number of bins for azimuthal averages", type=int, default=50
+)
 parser.add_argument(
     "-x", help="x position of the central point", type=float, default=1.0
 )
@@ -66,10 +77,7 @@ parser.add_argument(
 parser.add_argument(
     "-z", help="z position of the central point", type=float, default=1.0
 )
-parser.add_argument(
+
-    "--fft", help="use quick and dirty fft rendering", action="store_true"
-)
-parser.add_argument("-i", "--interactive", help="Interactive mode", action="store_true")
 
 parser.add_argument(
     "--colormap", help="Colormap used", choices=dp.P.colormaps(), default="plasma"
@@ -131,7 +139,7 @@ for run in runs:
                         fft=args.fft,
                         interactive=args.interactive,
                     )
-                if args.disk:
+                if args.disk or args.compare or args.evolution:
                     dp.disk_prop(
                         path_in,
                         i,
@@ -141,6 +149,7 @@ for run in runs:
                         force=args.force_redo,
                         pos_star=np.array([args.x, args.y, args.z]),
                     )
+                if args.disk:
                     dp.plot_disk_prop(
                         path_out,
                         i,
@@ -153,7 +162,9 @@ for run in runs:
                 if args.watch and failures < args.allowed_failures:
                     failures = failures + 1
                     print(
-                        "Unable to proceed for run {} output {}. Trying again in {} s ({} tries remaining)".format(
+                        "Unable to proceed for run {} \
+                    output {}. Trying again in {} s ({} \
+                    tries remaining)".format(
                             run, i, args.waiting_time, args.allowed_failures - failures
                         )
                     )
@@ -162,3 +173,17 @@ for run in runs:
                     break
                 else:
                     raise
+    if args.evolution:
+        dp.evolution(
+            path_out,
+            range(first, last + 1, step),
+            force=args.force_redo,
+            interactive=args.interactive,
+        )
+
+if args.compare:
+    path_suffix = project
+    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)