some updates to the massive star slice plots (#2759)

now we offset t based on the driving time
also add a time-sequence plot

Exec/science/massive_star/analysis/massive_star_multi.py

@@ -22,6 +22,13 @@ def make_plot(plotfile, prefix="plot", width_frac=0.1,
 
     ds = CastroDataset(plotfile)
 
+    t_drive = 0.0
+    if "[*] castro.drive_initial_convection_tmax" in ds.parameters:
+        t_drive = ds.parameters["[*] castro.drive_initial_convection_tmax"]
+    elif "castro.drive_initial_convection_tmax" in ds.parameters:
+        t_drive = ds.parameters["castro.drive_initial_convection_tmax"]
+    print(t_drive)
+
     xmin = ds.domain_left_edge[0]
     xmax = ds.domain_right_edge[0]
 
@@ -70,32 +77,33 @@ def make_plot(plotfile, prefix="plot", width_frac=0.1,
 
         if f == "enuc":
             # now do a contour of NSE
-            sp.annotate_contour("in_nse", ncont=1, clim=(0.5, 0.5), take_log=False,
-                                plot_args={"colors": "k", "linewidths": 2})
+            if ("boxlib", "in_nse") in ds.derived_field_list:
+                sp.annotate_contour("in_nse", ncont=1, clim=(0.5, 0.5), take_log=False,
+                                    plot_args={"colors": "k", "linewidths": 2})
 
     sp.set_axes_unit("cm")
 
     fig = sp.export_to_mpl_figure((layout[0], layout[1]), axes_pad=(1.0, 0.4),
                                   cbar_location=cbar_location, cbar_pad="2%")
 
     fig.subplots_adjust(left=0.05, right=0.95, bottom=0.025, top=0.975)
-    fig.text(0.02, 0.02, f"time = {float(ds.current_time):8.2f} s",
+    fig.text(0.02, 0.02, f"$t - \\tau_\\mathrm{{drive}}$ = {float(ds.current_time) - t_drive:6.1f} s",
              transform=fig.transFigure)
     fig.set_size_inches(size)
     fig.tight_layout()
     extra = ""
     if layout[0] >= layout[1]:
         extra = "_vertical"
 
-    fig.savefig(f"{prefix}_{os.path.basename(plotfile)}_w{width_frac:04.2f}{extra}.pdf", pad_inches=0)
+    fig.savefig(f"{prefix}_{os.path.basename(plotfile)}_w{width_frac:04.2f}{extra}.pdf", pad_inches=0.1, bbox_inches="tight")
 
 
 if __name__ == "__main__":
 
     p = argparse.ArgumentParser()
     p.add_argument("--vertical", action="store_true",
                    help="plot 2x2 or 1x4")
-    p.add_argument("--width_fraction", type=float, default=0.1,
+    p.add_argument("--width-fraction", type=float, default=0.1,
                    help="fraction of domain to show")
     p.add_argument("plotfile", type=str, nargs=1,
                    help="plotfile to plot")

Exec/science/massive_star/analysis/massive_star_sequence.py

@@ -0,0 +1,178 @@
+#!/usr/bin/env python3
+
+import argparse
+import os
+
+import matplotlib
+import matplotlib.pyplot as plt
+import numpy as np
+from mpl_toolkits.axes_grid1 import ImageGrid
+
+import yt
+from yt.frontends.boxlib.api import CastroDataset
+# assume that our data is in CGS
+from yt.units import cm
+
+matplotlib.use('agg')
+
+times = [50, 100, 150, 200, 250]
+
+def find_files(plist):
+
+    mask = np.zeros(len(times))
+    files_to_plot = []
+    for pfile in plist:
+        for k, t in enumerate(times):
+            if mask[k]:
+                continue
+            print(pfile)
+            ds = CastroDataset(pfile)
+            if ds.current_time >= t:
+                files_to_plot.append(pfile)
+                mask[k] = 1.0
+
+    return files_to_plot
+
+
+def make_plot(field, pfiles, *,
+              width_frac=0.1,
+              size=(19.2, 10.8)):
+
+    fig = plt.figure()
+
+    if len(pfiles) > 8:
+        nrows = 3
+        ncols = (len(pfiles) + 1)//3
+    elif len(pfiles) >= 4:
+        nrows = 2
+        ncols = (len(pfiles) + 1)//2
+    else:
+        nrows = 1
+        ncols = len(pfiles)
+
+    grid = ImageGrid(fig, 111, nrows_ncols=(nrows, ncols),
+                     axes_pad=0.2, cbar_pad=0.05,
+                     label_mode="L", cbar_mode="single")
+
+    # get some metadata
+    ds = CastroDataset(pfiles[0])
+
+    t_drive = 0.0
+    if "[*] castro.drive_initial_convection_tmax" in ds.parameters:
+        t_drive = ds.parameters["[*] castro.drive_initial_convection_tmax"]
+    elif "castro.drive_initial_convection_tmax" in ds.parameters:
+        t_drive = ds.parameters["castro.drive_initial_convection_tmax"]
+    print(t_drive)
+
+    xmin = ds.domain_left_edge[0]
+    xmax = ds.domain_right_edge[0]
+
+    ymin = ds.domain_left_edge[1]
+    ymax = ds.domain_right_edge[1]
+
+    xctr = 0.0 * xmin
+    L_x = xmax - xmin
+
+    yctr = 0.5 * (ymin + ymax)
+    L_y = ymax - ymin
+
+    f = field
+
+    if f == "Temp" or f == "abar":
+        text_color = "white"
+    else:
+        text_color = "black"
+
+    for i in range(nrows * ncols):
+
+        if i < len(pfiles):
+            pf = pfiles[i]
+        else:
+            grid[i].remove()
+            continue
+
+        ds = CastroDataset(pf)
+
+        sp = yt.SlicePlot(ds, "theta", f,
+                          center=[xmin + 0.5*width_frac*L_x, yctr, 0.0*cm],
+                          width=[width_frac*L_x, width_frac*L_y, 0.0*cm],
+                          fontsize="12")
+
+        sp.set_buff_size((2400,2400))
+
+        sp.annotate_text((0.05, 0.05),
+                         f"$t - \\tau_\\mathrm{{drive}}$ = {float(ds.current_time) - t_drive:8.2f} s",
+                         coord_system="axis",
+                         text_args={"color": text_color})
+
+
+        if f == "Ye":
+            sp.set_zlim(f, 0.46, 0.5)
+            sp.set_log(f, False)
+            sp.set_cmap(f, "magma_r")
+        elif f == "abar":
+            sp.set_log(f, False)
+            sp.set_cmap(f, "viridis")
+        elif f == "enuc":
+            sp.set_log(f, True, linthresh=1.e12)
+            sp.set_zlim(f, -1.e20, 1.e20)
+            sp.set_cmap(f, "bwr")
+        elif f == "MachNumber":
+            sp.set_zlim(f, 1.e-4, 0.3)
+            sp.set_cmap(f, "plasma")
+        elif f == "magvel":
+            sp.set_zlim(f, 100.0, 2.e7)
+            sp.set_cmap(f, "viridis")
+        elif f == "magvort":
+            sp.set_cmap(f, "magma")
+            sp.set_zlim(f, 1.e-2, 5)
+
+        if f == "enuc":
+            # now do a contour of NSE
+            if ("boxlib", "in_nse") in ds.derived_field_list:
+                sp.annotate_contour("in_nse", levels=1, clim=(0.5, 0.5), take_log=False,
+                                    plot_args={"colors": "k", "linewidths": 2})
+
+        sp.set_axes_unit("cm")
+
+        plot = sp.plots[field]
+        plot.figure = fig
+        plot.axes = grid[i].axes
+        plot.cax = grid.cbar_axes[i]
+        if i < len(pfiles)-1:
+            grid[i].axes.xaxis.offsetText.set_visible(False)
+        if i > 0:
+            grid[i].axes.yaxis.offsetText.set_visible(False)
+
+        sp._setup_plots()
+
+    fig.set_size_inches(size)
+    fig.tight_layout()
+    extra = ""
+
+    fig.savefig(f"massive_star_{f}_w{width_frac:04.2f}.pdf", pad_inches=0.1, bbox_inches="tight")
+
+
+if __name__ == "__main__":
+
+    p = argparse.ArgumentParser()
+    p.add_argument("--var", type=str, default="Temp",
+                   help="variable to plot")
+    p.add_argument("--vertical", action="store_true",
+                   help="plot 2x2 or 1x4")
+    p.add_argument("--width-fraction", type=float, default=0.1,
+                   help="fraction of domain to show")
+    p.add_argument("plotfiles", type=str, nargs="+",
+                   help="list of plotfiles to plot")
+
+    args = p.parse_args()
+
+    if args.vertical:
+        size = (7.5, 8.5)
+    else:
+        size = (19.2, 8.5)
+
+    plist = find_files(args.plotfiles)
+
+    make_plot(args.var, plist,
+              width_frac=args.width_fraction, size=size)