improvements to xrb plotting script (#3154)

this is for xrb_spherical

Author: zhichen3

Exec/science/xrb_spherical/GNUmakefile

@@ -26,7 +26,7 @@ NUM_MODELS := 2
 EOS_DIR     := helmholtz
 
 # This sets the network directory in $(MICROPHYSICS_HOME)/networks
-NETWORK_DIR := he-burn/he-burn-18a
+NETWORK_DIR := he-burn/he-burn-19am
 
 INTEGRATOR_DIR := RKC
 

Exec/science/xrb_spherical/Problem_Derive.H

@@ -0,0 +1 @@
+../flame_wave/Problem_Derive.H

Exec/science/xrb_spherical/Problem_Derive.cpp

@@ -0,0 +1 @@
+../flame_wave/Problem_Derive.cpp

Exec/science/xrb_spherical/Problem_Derives.H

@@ -0,0 +1 @@
+../flame_wave/Problem_Derives.H

Exec/science/xrb_spherical/analysis/flame_speed.py

@@ -28,7 +28,7 @@
 from front_tracker.py to plot the time evolution of flame front θ.
 ''')
 
-parser.add_argument('tracking_fname', type=str,
+parser.add_argument('tracking_fnames', nargs='+', type=str,
                     help='cvs file generated from front_tracker.py to track flame front position')
 parser.add_argument('--tmin', default=0.0, type=float,
                     help='minimum time for curve fitting. Note that this will not affect plotting')
@@ -37,20 +37,29 @@
 
 args = parser.parse_args()
 
-# Read in data
-# data has columns: fname, time, front_theta, theta_max_avg, max_avg, theta_max, max_val.
-tracking_data = pd.read_csv(args.tracking_fname)
+all_data = []
+
+# Loop over all tracking files and append them
+for fname in args.tracking_fnames:
+    df = pd.read_csv(fname)
+    all_data.append(df)
+
+# concatenate all files together
+tracking_data = pd.concat(all_data, ignore_index=True)
 
+# sort by the time
+tracking_data = tracking_data.sort_values(by='time[ms]')
+
+# data has columns: fname, time, front_theta, theta_max_avg, max_avg, theta_max, max_val.
 # Get time and theta, these should already be time-sorted
 times = tracking_data['time[ms]']
 front_thetas = tracking_data['front_theta']
 
 # Only plot up to tmax
 if args.tmax is not None:
     cond = times < args.tmax
-
-times = times[cond]
-front_thetas = front_thetas[cond]
+    times = times[cond]
+    front_thetas = front_thetas[cond]
 
 # Now do a curve fit to the data.
 # Now apply tmin so that we ignore the transient phase during fitting

Exec/science/xrb_spherical/analysis/inset_slice.py

@@ -19,9 +19,8 @@
 def single_slice(ds, field:str,
                  loc: str = "top", widthScale: float = 3.0,
                  theta: Optional[float] = None,
+                 position: float | None = None,
                  displace_theta: bool = True,
-                 annotate_vline: bool = True,
-                 annotate_lat_lines: bool = True,
                  show_full_star: bool = True) -> None:
     """
     A slice plot a single dataset for a single field parameters for Spherical2D geometry.
@@ -45,14 +44,13 @@ def single_slice(ds, field:str,
     theta:
       user defined theta center of the slice plot
 
+    position:
+      draw a vertical line on user defined front position in theta (in radian)
+
     displace_theta:
       When theta is explicitly defined, do we want to displace the center
       of the slice plot by ~0.7. This is helpful when tracking the flame front.
 
-    annotate_vline:
-      do we want to annotate a vertical line at where theta is.
-      This is used to indicate flame front.
-
     show_full_star:
       do we want to plot the full star instead of a zoom-in slice plot.
     """
@@ -96,18 +94,17 @@ def single_slice(ds, field:str,
     # sp.annotate_text((0.05, 0.05), f"{currentTime.in_cgs():8.5f} s")
 
     # Plot a vertical to indicate flame front
-    if theta is not None and annotate_vline:
-        sp.annotate_line([r[0]*np.sin(theta), r[0]*np.cos(theta)],
-                         [r[2]*np.sin(theta), r[2]*np.cos(theta)],
+    if position is not None:
+        sp.annotate_line([r[0]*np.sin(position), r[0]*np.cos(position)],
+                         [r[2]*np.sin(position), r[2]*np.cos(position)],
                          coord_system="plot",
                          color="k",
                          linewidth=1.5,
                          linestyle="-.")
 
     ### Annotate Latitude Lines
-    if annotate_lat_lines:
-        annotate_latitude_lines(sp, center, box_widths, r,
-                                show_full_star=show_full_star)
+    annotate_latitude_lines(sp, center, box_widths, r,
+                            show_full_star=show_full_star)
 
     sp._setup_plots()
     return sp
@@ -131,13 +128,14 @@ def single_slice(ds, field:str,
     parser.add_argument('-t', '--theta', type=float,
                         help="""user defined theta center location of the plot domain.
                         Alternative way of defining plotting center""")
+    parser.add_argument('-p', '--position', type=float,
+                        help="""user defined front position in theta,
+                        this will draw a vertical line to annotate the front position""")
     parser.add_argument('-w', '--width', default=2.0, type=float,
                         help="scaling for the domain width of the slice plot")
     parser.add_argument('--displace_theta', action='store_true',
                         help="""whether to displace the theta that defines the center of the frame.
                         This is useful when theta represents the flame front position.""")
-    parser.add_argument('--annotate_vline', action='store_true',
-                        help="whether to annotate a vertical line along the given theta")
 
     args = parser.parse_args()
 
@@ -152,9 +150,9 @@ def single_slice(ds, field:str,
 
     # First get the slice plot of the full-star.
     full_star_slice = single_slice(ds, args.field, loc=loc,
-                                   widthScale=args.width, theta=args.theta,
-                                   displace_theta=args.displace_theta, annotate_vline=args.annotate_vline,
-                                   annotate_lat_lines=True, show_full_star=True)
+                                   widthScale=args.width, theta=args.theta, position=args.position,
+                                   displace_theta=args.displace_theta,
+                                   show_full_star=True)
     # full_star_slice.render()
 
     # Extract the figure of the full-star slice and use that as the main figure for plotting.
@@ -166,9 +164,9 @@ def single_slice(ds, field:str,
 
     # Get the slice of the zoom-in plot
     zoom_in_slice = single_slice(ds, args.field, loc=loc,
-                                 widthScale=args.width, theta=args.theta,
-                                 displace_theta=args.displace_theta, annotate_vline=args.annotate_vline,
-                                 annotate_lat_lines=True, show_full_star=False)
+                                 widthScale=args.width, theta=args.theta, position=args.position,
+                                 displace_theta=args.displace_theta,
+                                 show_full_star=False)
     zoom_in_slice.hide_colorbar()
     # zoom_in_slice.render()
 
@@ -226,10 +224,21 @@ def single_slice(ds, field:str,
 
     ### Now annotate inset box lines ###
     # loc1 and loc2: corners to connect (1: upper right, 2: upper left, 3: lower left, 4: lower right)
-    if args.theta < np.pi/3.0 or (args.theta is None and loc=="top"):
+    if args.theta is None:
+        if loc == "top":
+            loc1 = 1
+            loc2 = 2
+        elif loc == "bot":
+            loc1 = 3
+            loc2 = 4
+        else:
+            # mid
+            loc1 = 1
+            loc2 = 4
+    elif args.theta < np.pi/3.0:
         loc1 = 1
         loc2 = 2
-    elif args.theta > 2.0*np.pi/3.0 or (args.theta is None and loc=="bot"):
+    elif args.theta > 2.0*np.pi/3.0:
         loc1 = 3
         loc2 = 4
     else:

Exec/science/xrb_spherical/analysis/slice.py

@@ -7,11 +7,41 @@
 import math
 from typing import List, Optional
 import numpy as np
+import pynucastro as pyna
 import matplotlib.pyplot as plt
 from mpl_toolkits.axes_grid1 import ImageGrid
 from yt.frontends.boxlib.api import CastroDataset
 from yt.units import km
 
+def _ash(field, data):
+    '''
+    Computes the rho X_ash.
+    Here ash is anything heavier than Oxygen, but exclude Fe and Ni
+    '''
+
+    ds = data.ds
+    field_list = ds.field_list
+
+    # If X(ash) is directly stored, then just use that
+    if ("boxlib", "X(ash)") in field_list:
+        return data["boxlib", "X(ash)"] * data["boxlib", "density"]
+
+    # If we cannot find X(ash) as a available  field then compute manually
+    rhoAsh = 0
+    for f in field_list:
+        # If the first two letters are "X(", then we're dealing with species massfractions
+        if f[1][:2] == "X(":
+            # Then extract out what species we have, assuming the format is "X(...)"
+            speciesName = f[1][2:-1]
+            nuc = pyna.Nucleus(speciesName)
+
+            # Include elements beyond oxygen but don't include Ni56
+            if nuc.Z > 8.0 and nuc.Z != 56:
+                rhoAsh += data["boxlib", "density"] * data[f]
+
+    return rhoAsh
+
+
 def extract_info(ds,
                  loc: str = "top", widthScale: float = 3.0,
                  widthRatio: float = 1.0,
@@ -220,9 +250,8 @@ def slice(fnames:list[str], fields:list[str],
           loc: str = "top", widthScale: float = 3.0,
           widthRatio: float = 1.0,
           theta: float | None = None,
+          position: float | None = None,
           displace_theta: bool = False,
-          annotate_vline: bool = False,
-          annotate_lat_lines: bool = True,
           show_full_star: bool = False) -> None:
     """
     A slice plot of the datasets for different field parameters for Spherical2D geometry.
@@ -252,19 +281,15 @@ def slice(fnames:list[str], fields:list[str],
       For widthRatio > 1, the vertical width is larger than horizontal.
 
     theta:
-      user defined theta center of the slice plot
+      user defined theta (in radian) center of the slice plot
+
+    position:
+      draw a vertical line on user defined front position in theta (in radian)
 
     displace_theta:
       whether to displace theta so that the vertical lines that represents
       the flame front is offset by some amount
 
-    annotate_vline:
-      whether to plot a vertical line to represent the flame front,
-      which is represented by what theta is.
-
-    annotate_lat_lines:
-      whether to annotate latitude lines.
-
     show_full_star:
       whether to plot the full star rather than a zoom-in
     """
@@ -278,7 +303,7 @@ def slice(fnames:list[str], fields:list[str],
     nx = math.ceil(num/ny)
 
     grid = ImageGrid(fig, 111, nrows_ncols=(nx, ny),
-                     axes_pad=1, label_mode="L", cbar_location="right",
+                     axes_pad=1.1, label_mode="L", cbar_location="right",
                      cbar_mode="each", cbar_size="2.5%", cbar_pad="0%")
 
     # Output plot file name
@@ -292,6 +317,12 @@ def slice(fnames:list[str], fields:list[str],
                                              theta=theta,
                                              displace_theta=displace_theta,
                                              show_full_star=show_full_star)
+
+        #add rhoX_ash as a derived field
+        ds.add_field(("gas", "ash"), function=_ash,
+                     display_name=r"\rho X\left(ash\right)",
+                     units="auto", sampling_type="cell")
+
         for i, field in enumerate(fields):
             # Plot each field parameter
             sp = yt.SlicePlot(ds, 'phi', field, width=box_widths, fontsize=20)
@@ -310,6 +341,10 @@ def slice(fnames:list[str], fields:list[str],
                 sp.set_zlim(field, 4, 8)
                 sp.set_log(field, False)
                 sp.set_cmap(field, "plasma_r")
+            elif field == "ash":
+                sp.set_zlim(field, 1.e-2, 1e6)
+                sp.set_log(field, True)
+                sp.set_cmap(field, "plasma_r")
             elif field == "enuc":
                 sp.set_zlim(field, 1.e15, 1.e20)
                 sp.set_log(field, linthresh=1.e11)
@@ -321,18 +356,17 @@ def slice(fnames:list[str], fields:list[str],
             # sp.annotate_text((0.05, 0.05), f"{currentTime.in_cgs():8.5f} s")
 
             # Plot a vertical to indicate flame front
-            if theta is not None and annotate_vline:
-                sp.annotate_line([r[0]*np.sin(theta), r[0]*np.cos(theta)],
-                                 [r[2]*np.sin(theta), r[2]*np.cos(theta)],
+            if position is not None:
+                sp.annotate_line([r[0]*np.sin(position), r[0]*np.cos(position)],
+                                 [r[2]*np.sin(position), r[2]*np.cos(position)],
                                  coord_system="plot",
                                  color="k",
                                  linewidth=1.5,
                                  linestyle="-.")
 
             ### Annotate Latitude Lines
-            if annotate_lat_lines:
-                annotate_latitude_lines(sp, center, box_widths, r,
-                                        show_full_star=show_full_star)
+            annotate_latitude_lines(sp, center, box_widths, r,
+                                    show_full_star=show_full_star)
 
             plot = sp.plots[field]
             plot.figure = fig
@@ -391,6 +425,9 @@ def slice(fnames:list[str], fields:list[str],
     parser.add_argument('-t', '--theta', type=float,
                         help="""user defined theta center location of the plot domain.
                         Alternative way of defining plotting center""")
+    parser.add_argument('-p', '--position', type=float,
+                        help="""user defined front position in theta,
+                        this will draw a vertical line to annotate the front position""")
     parser.add_argument('-r', '--ratio', default=1.0, type=float,
                         help="""The ratio between the horizontal and vertical width of the slice plot.
                         For ratio < 1, horizontal width is larger than vertical.
@@ -400,8 +437,6 @@ def slice(fnames:list[str], fields:list[str],
     parser.add_argument('--displace_theta', action='store_true',
                         help="""whether to displace the theta that defines the center of the frame.
                         This is useful when theta represents the flame front position.""")
-    parser.add_argument('--annotate_vline', action='store_true',
-                        help="whether to annotate a vertical line along the given theta")
     parser.add_argument('--show_full_star', action='store_true',
                         help="whether show the full star in the background")
 
@@ -417,6 +452,7 @@ def slice(fnames:list[str], fields:list[str],
         parser.error("loc must be one of the three: {top, mid, bot}")
 
     slice(args.fnames, args.fields, loc=loc,
-          widthScale=args.width, widthRatio=args.ratio, theta=args.theta,
-          displace_theta=args.displace_theta, annotate_vline=args.annotate_vline,
-          annotate_lat_lines=True, show_full_star=args.show_full_star)
+          widthScale=args.width, widthRatio=args.ratio,
+          theta=args.theta, position=args.position,
+          displace_theta=args.displace_theta,
+          show_full_star=args.show_full_star)

Exec/science/xrb_spherical/analysis/slice_sequence.py

@@ -25,8 +25,6 @@
 parser.add_argument('--displace_theta', action='store_true',
                     help="""whether to displace the theta that defines the center of the frame.
                         This is useful when theta represents the flame front position.""")
-parser.add_argument('--annotate_vline', action='store_true',
-                    help="whether to annotate a vertical line along the given theta")
 parser.add_argument('--jobs', '-j', default=1, type=int,
                     help="""Number of workers to plot in parallel""")
 
@@ -40,13 +38,16 @@
 fnames = tracking_data["fname"]
 front_thetas = tracking_data["front_theta"]
 
+# Process front_thetas with moving average so that frame doesn't jitter
+window = 10
+smooth_thetas = tracking_data['front_theta'].rolling(window, center=True).mean()
+
 # Parallelize the plotting
 with ProcessPoolExecutor(max_workers=args.jobs) as executor:
     future_to_index = {
         executor.submit(slice, [fname], args.fields, widthScale=args.width,
-                        widthRatio=args.ratio, theta=front_thetas[i],
-                        displace_theta=args.displace_theta, annotate_vline=args.annotate_vline,
-                        annotate_lat_lines=True): i
+                        widthRatio=args.ratio, theta=smooth_thetas[i], position=front_thetas[i],
+                        displace_theta=args.displace_theta): i
         for i, fname in enumerate(fnames)
     }
     try:

…nce/xrb_spherical/inputs.He.10deg.1000Hz …_spherical/inputs/inputs.He.10deg.1000HzExec/science/xrb_spherical/inputs.He.10deg.1000Hz renamed to Exec/science/xrb_spherical/inputs/inputs.He.10deg.1000Hz Exec/science/xrb_spherical/inputs.He.10deg.1000Hz renamed to Exec/science/xrb_spherical/inputs/inputs.He.10deg.1000Hz

@@ -107,7 +107,9 @@ amr.derive_plot_vars = ALL
 amr.small_plot_file         = xrb_spherical_1000Hz_smallplt        # root name of plotfile
 amr.small_plot_per          = 1.e-4      # number of seconds between plotfiles
 amr.small_plot_vars = density Temp
-amr.derive_small_plot_vars = abar x_velocity y_velocity z_velocity enuc
+amr.derive_small_plot_vars = abar x_velocity y_velocity z_velocity enuc X(ash)
+
+amr.file_name_digits = 7
 
 # problem initialization
 

…ce/xrb_spherical/inputs.He.180deg.1000Hz …spherical/inputs/inputs.He.180deg.1000HzExec/science/xrb_spherical/inputs.He.180deg.1000Hz renamed to Exec/science/xrb_spherical/inputs/inputs.He.180deg.1000Hz Exec/science/xrb_spherical/inputs.He.180deg.1000Hz renamed to Exec/science/xrb_spherical/inputs/inputs.He.180deg.1000Hz

@@ -107,7 +107,9 @@ amr.derive_plot_vars = ALL
 amr.small_plot_file         = xrb_spherical_1000Hz_smallplt        # root name of plotfile
 amr.small_plot_per          = 1.e-4      # number of seconds between plotfiles
 amr.small_plot_vars = density Temp
-amr.derive_small_plot_vars = abar x_velocity y_velocity z_velocity enuc
+amr.derive_small_plot_vars = abar x_velocity y_velocity z_velocity enuc X(ash)
+
+amr.file_name_digits = 7
 
 # problem initialization