Add map for Q

disk_postprocess.py

@@ -19,8 +19,11 @@ from pymses.analysis import Camera, raytracing, slicing
 from pymses.filters import CellsToPoints
 from pymses.analysis import ScalarOperator, FractionOperator, MaxLevelOperator
 
+# extension for out files
+out_ext = ".jpeg"
+
 P.rcParams["image.cmap"] = "plasma"
-P.rcParams["savefig.dpi"] = 800
+P.rcParams["savefig.dpi"] = 400
 
 
 def make_image_disk(
@@ -52,7 +55,7 @@ def make_image_disk(
     """
 
     ro = pymses.RamsesOutput(path, num, order=order)
-    amr = ro.amr_source(["rho", "vel", "P", "Bl", "Br"])
+    amr = ro.amr_source(["rho", "vel", "P"])
     rad = 0.5
     center = [0.5, 0.5, 0.5]
 
@@ -86,6 +89,7 @@ def make_image_aux(
     vel_red=20,
     tag="",
     cpuamr=False,
+    pos_star=np.array([1.0, 1.0, 1.0]),
     put_title=True,
 ):
     """
@@ -122,7 +126,7 @@ def make_image_aux(
     else:
         directory = path
 
-    name = directory + "/coldens_z" + "_" + tag + "_" + format(num, "05") + ".jpeg"
+    name = directory + "/coldens_z" + "_" + tag + "_" + format(num, "05") + out_ext
     if len(glob.glob(name)) == 1 and not force:
         return
 
@@ -181,7 +185,7 @@ def make_image_aux(
         cbar = P.colorbar(im)
         cbar.set_label(r"$log(N)$ (code)")
         name = directory + "/coldens_" + ax_los + "_" + tag + "_" + format(num, "05")
-        name_im = name + ".jpeg"
+        name_im = name + out_ext
         P.savefig(name_im)
 
         # Rho slice
@@ -248,7 +252,7 @@ def make_image_aux(
         cbar.set_label(r"$log(n)$ (code)")
 
         name = directory + "/rho_" + ax_los + "_" + tag + "_" + format(num, "05")
-        name_im = name + ".jpeg"
+        name_im = name + out_ext
         P.savefig(name_im)
 
         P_op = ScalarOperator(lambda dset: dset["P"], ro.info["unit_pressure"])
@@ -279,18 +283,48 @@ def make_image_aux(
         cbar.set_label(r"$log(T) \, (K)$")
 
         name = directory + "/T_" + ax_los + "_" + tag + "_" + format(num, "05")
-        name_im = name + ".jpeg"
+        name_im = name + out_ext
         P.savefig(name_im)
 
         # Toomre parameter
 
-        if False and ax_los == "z":
+        if ax_los == "z":
-            map_Q = np.log10(np.sqrt(dmap_P.map.T / dmap_rho.map.T)) / (
-                np.pi * map_col * G
-            )
 
+            def omega_rho_func(dset):
+                pos = dset.get_cell_centers()
+                pos = pos - (pos_star / lbox)
+
+                xx = pos[:, :, 0]
+                yy = pos[:, :, 1]
+                rc = np.sqrt(xx ** 2 + yy ** 2)  # cylindrical radius
+                vx = dset["vel"][:, :, 0]
+                vy = dset["vel"][:, :, 1]
+                omega_rho = 1.0 / rc ** 2
+                omega_rho = omega_rho * dset["rho"]
+                vyx = vy * xx
+                vxy = vx * yy
+                omega_rho = omega_rho * (vyx - vxy)
+                return omega_rho
+
+            omega_op = FractionOperator(
+                omega_rho_func, lambda dset: dset["rho"], 1.0 / ro.info["unit_time"]
+            )
+            cs_op = FractionOperator(
+                lambda dset: dset["P"],
+                lambda dset: dset["rho"],
+                ro.info["unit_velocity"],
+            )
+            rt_omega = raytracing.RayTracer(amr, ro.info, omega_op)
+            rt_cs = raytracing.RayTracer(amr, ro.info, cs_op)
+
+            dmap_omega = rt_omega.process(cam)
+            dmap_cs = rt_cs.process(cam)
+            dmap_Q = (lbox * dmap_cs.map.T) * dmap_omega.map.T / (np.pi * G * dmap_col)
+            map_Q = np.log10(dmap_Q)
+
+            P.close()
             im = P.imshow(
-                map_T,
+                map_Q,
                 extent=[
                     (-radius + center[0]) * lbox_units,
                     (radius + center[0]) * lbox_units,
@@ -305,14 +339,14 @@ def make_image_aux(
 
             if put_title:
                 P.title(title)
-                P.xlabel(title_ax[ax_h])
+            P.xlabel(title_ax[ax_h])
-                P.ylabel(title_ax[ax_v])
+            P.ylabel(title_ax[ax_v])
-                cbar = P.colorbar(im)
+            cbar = P.colorbar(im)
-                cbar.set_label(r"$log(T) \, (K)$")
+            cbar.set_label(r"$log(Q)$")
 
-                name = directory + "/Q_" + ax_los + "_" + tag + "_" + format(num, "05")
+            name = directory + "/Q_" + ax_los + "_" + tag + "_" + format(num, "05")
-                name_im = name + ".jpeg"
+            name_im = name + out_ext
-                P.savefig(name_im)
+            P.savefig(name_im)
 
         if cpuamr:
             level_op = MaxLevelOperator()
@@ -344,7 +378,7 @@ def make_image_aux(
             cbar.set_label(r"level")
 
             name = directory + "/level_" + ax_los + "_" + tag + "_" + format(num, "05")
-            name_im = name + ".jpeg"
+            name_im = name + out_ext
             P.savefig(name_im)
 
             cpu_op = ScalarOperator(
@@ -378,7 +412,7 @@ def make_image_aux(
             cbar.set_label(r"cpu")
 
             name = directory + "/cpu_" + ax_los + "_" + tag + "_" + format(num, "05")
-            name_im = name + ".jpeg"
+            name_im = name + out_ext
             P.savefig(name_im)
 
 
@@ -432,7 +466,7 @@ def disk_prop(
     ro = pymses.RamsesOutput(path_in, num)
     lbox = ro.info["boxlen"]  # boxlen in codeunits (=>pc)
 
-    time = ro.info["time"]  # * scale_t / Myr
+    time = ro.info["time"]  # time in codeunits
 
     # Get array of cell positions
     amr = ro.amr_source(["rho", "vel", "Br", "Bl", "P"])
@@ -440,9 +474,9 @@ def disk_prop(
     cells = cell_source.flatten()
     dx = cells.get_sizes()
     pos = cells.points * lbox
-
     # Get positions in the frame of the protostar
     pos = pos - pos_star
+
     # Get cylindrical radius
     rc = np.sqrt(pos[:, 0] ** 2 + pos[:, 1] ** 2)
     # Get velocities
@@ -474,6 +508,8 @@ def disk_prop(
     v_az_disk = v_az[mask]
     v_kepl = np.sqrt(mass_star * G / rc_disk)
 
+    total_mass_disk = np.sum(rho_disk * dvol_disk)
+
     # Initialize binned quantities
     cs_rad = np.zeros(nb_bin - 1)
     temp_rad = np.zeros(nb_bin - 1)
@@ -552,6 +588,7 @@ def disk_prop(
 
     prop_disk = {
         "time": time,
+        "tot_mass": total_mass_disk,
         "rad": rad[0 : nb_bin - 1],
         "center": pos_star,
         "alpha_rey": alpha_rey_rad,
@@ -572,16 +609,13 @@ def disk_prop(
     f.close()
 
 
-def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
+def plot_disk_prop(path, num, force=False, tag=""):
     """
     Plot properties of a disk
 
     num       id of the ramses output
     path      path to the properties file
     force     if set, redo plots even if already done
-
-    pdf       if set, do output in pdf as well
-
     """
 
     # Load property file
@@ -595,16 +629,22 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
     f.close()
 
     # Check if the output file exists, and exit if it is the case
-    name_save = path + "/rho_disk_r_" + str(num).zfill(5) + ".jpeg"
+    name_save = path + "/rho_disk_r_" + str(num).zfill(5) + out_ext
     if not force and len(glob.glob(name_save)) != 0:
         return
 
+    time = prop_disk["time"]
+    mass = prop_disk["tot_mass"]
+    title = "t=" + str(time)[0:5] + " (code)"
+
     P.close()
     P.plot(
         np.log10(prop_disk["rad"]), np.log10(prop_disk["rho"]), color="k", linewidth=2
     )
-    P.ylabel(r"$\log(n) \, (cm^{-3})$")
+    P.ylabel(r"$\log(n) \, (code)$")
     P.xlabel("disk radius")
+    P.title(title)
+    P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + out_ext)
 
     if pdf:
         P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + ".pdf")
@@ -616,6 +656,8 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
     )
     P.ylabel(r"$\log(T) \, (K)$")
     P.xlabel("disk radius")
+    P.title(title)
+    P.savefig(path + "/T_disk_r_" + str(num).zfill(5) + out_ext)
 
     if pdf:
         P.savefig(path + "/T_disk_r_" + str(num).zfill(5) + ".pdf")
@@ -636,9 +678,7 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
     P.ylabel(r"$V \, (km s^{-1})$")
     P.xlabel("disk radius")
 
-    if pdf:
+    P.savefig(path + "/V_disk_r_" + str(num).zfill(5) + out_ext)
-        P.savefig(path + "/V_disk_r_" + str(num).zfill(5) + ".pdf")
-    P.savefig(path + "/V_disk_r_" + str(num).zfill(5) + ".jpeg")
 
     P.close()
     P.plot(
@@ -649,10 +689,8 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
     )
     P.ylabel(r"$\log(N) \, (cm^{-2})$")
     P.xlabel("disk radius ")
-
+    P.title(title)
-    if pdf:
+    P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + out_ext)
-        P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + ".pdf")
-    P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + ".jpeg")
 
     # Alpha
     P.close()
@@ -664,10 +702,8 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
 
     P.ylabel(r"$\alpha$")
     P.xlabel("disk radius ")
-
+    P.title(title)
-    if pdf:
+    P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + out_ext)
-        P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + ".pdf")
-    P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + ".jpeg")
 
     # Q
 
@@ -675,10 +711,10 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
     P.xscale("log")
     P.yscale("log", linthreshy=0.001)
 
-    P.plot(prop_disk["rad"], abs(prop_disk["Q_kepl"]), color="b", linewidth=2)
-
     P.ylabel(r"$Q$")
     P.xlabel("disk radius ")
+    P.title(title + ", mass of disk = {} (code)".format(mass))
+    P.savefig(path + "/Q_r_" + str(num).zfill(5) + out_ext)
 
     if pdf:
         P.savefig(path + "/Q_r_" + str(num).zfill(5) + ".pdf")

pipeline_disk.py

@@ -14,18 +14,21 @@ storage_out = "/dsm/anais/storageA/"
 
 parser = argparse.ArgumentParser()
 parser.add_argument("runs", help="name of runs", nargs="*", default=["015_iso"])
-parser.add_argument(
+parser.add_argument("-b", "--begin", help="id of first output", type=int, default=1)
-    "-f", "--first_output", 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(
-    "-l", "--last_output", 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("-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"
 )
@@ -52,8 +55,8 @@ folder = "simus"
 project = args.project
 runs = args.runs
 
-first = args.first_output
+first = args.begin
-last = args.last_output
+last = args.end
 step = args.step
 
 
@@ -74,13 +77,23 @@ for run in runs:
         while not success:
             try:
                 dp.make_image_disk(
-                    path_in, i, path_out=path_out, tag=run, map_size=1024
+                    path_in,
+                    i,
+                    path_out=path_out,
+                    tag=run,
+                    map_size=1024,
+                    force=args.force_redo,
                 )
                 if args.disk:
                     dp.disk_prop(
-                        path_in, i, path_out=path_out, rad_ext=1, nb_bin=50, force=False
+                        path_in,
+                        i,
+                        path_out=path_out,
+                        rad_ext=1,
+                        nb_bin=50,
+                        force=args.force_redo,
                     )
-                    dp.plot_disk_prop(path_out, i, tag=run)
+                    dp.plot_disk_prop(path_out, i, tag=run, force=args.force_redo)
                 success = True
             except ValueError:
                 if args.watch and failures < args.allowed_failures: