Update xrb spherical problem

Exec/science/xrb_spherical/GNUmakefile

@@ -10,9 +10,9 @@ COMP       = gnu
 USE_MPI    = TRUE
 
 USE_GRAV   = TRUE
-USE_REACT  = FALSE
+USE_REACT  = TRUE
 
-USE_ROTATION  = FALSE
+USE_ROTATION  = TRUE
 USE_DIFFUSION = FALSE
 
 # define the location of the CASTRO top directory

Exec/science/xrb_spherical/analysis/front_tracker.py

@@ -0,0 +1,77 @@
+#!/usr/bin/env python3
+
+import sys
+import glob
+import yt
+import numpy as np
+from yt.frontends.boxlib.api import CastroDataset
+
+
+def track_flame_front(ds, metric):
+    '''
+    This function tracks the flame front position for a given dataset.
+    It returns a tuple of the form: (Time (in ms), Theta)
+
+    Procedure to determine flame front:
+    1) User selects e a quantity to use as a metric: enuc or Temp
+    2) Determine the global max of that quantity
+    3) Determine the theta position of the cell that contains the global max
+    4) Do a radial average of the data set to convert to 1D as a function of theta
+    5) Determine flame front at theta where the metric quantity drops to
+       metric_number * global_max of that quantity.
+    '''
+
+    time = ds.current_time.in_units("ms")
+
+
+    ad = ds.all_data()
+
+
+    # 1) Global max temperature: this can be used to track initial
+    # detonation resulted from the initial temperature perturbation
+    # 2) Use radially averaged enuc then find flame front is determined by
+    # when the averaged enuc drops below some percentage of max global enuc
+
+
+    return timeTheta
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="""
+    This file tracks the flame front and writes them into a txt file.
+    """)
+
+    parser.add_argument('--fnames', nargs='+', type=str,
+                        help="Dataset file names for tracking flame front.")
+    parser.add_argument('--field', '-f', default="enuc", type=str,
+                        metavar="{enuc, Temp}",
+                        help="""field parameter used as metric to determine flame front.
+                        Choose between {enuc, Temp}""")
+    parser.add_argument('--percent', '-p', default=0.001, type=float,
+                        help="""Float number between (0, 1]. Representing the percent of
+                        the global maximum of the field quantity used to track the flame.""")
+
+    args = parser.parse_args()
+
+    metric_quantities = ["enuc", "Temp"]
+    if args.field not in metric_quantities:
+        parser.error("field must be either enuc or Temp")
+
+    if args.percent <= 0.0 or args.percent > 1.0:
+        parser.error("metric must be a float between (0, 1]")
+
+    metric = (args.field, args.percent)
+    timeThetaArray = []
+
+    for fname in args.fnames:
+        ds = CastroDataset(fname)
+
+        # Get tuple in form (theta, time)
+        timeTheta= track_flame_front(ds, metric)
+        timeThetaArray.append(timeTheta)
+
+    # Sort array by time and write to file
+    timeThetaArray.sort()
+    timeThetaArray = np.array(timeThetaArray)
+
+    np.savetxt('front_tracking.dat', timeThetaArray, delimiter=',')

Exec/science/xrb_spherical/analysis/r_profile.py

@@ -41,6 +41,8 @@
     for theta in thetas:
         # simply go from (rmin, theta) -> (rmax, theta). Doesn't need to convert to physical R-Z
         ray = ds.ray((rmin, theta, 0*cm), (rmax, theta, 0*cm))
+
+        # sort by "t", which goes from 0 to 1 representing the spatial order.
         isrt = np.argsort(ray["t"])
         axes[i].plot(ray['r'][isrt], ray[f][isrt], label=r"$\theta$ = {:.4f}".format(float(theta)))
 

Exec/science/xrb_spherical/analysis/slice.py

@@ -3,92 +3,181 @@
 import sys
 import os
 import yt
+import argparse
+import math
+from typing import List, Optional
 import numpy as np
 import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import ImageGrid
 from yt.frontends.boxlib.api import CastroDataset
 
-from yt.units import cm
 
-"""
-Given a plot file and field name, it gives slice plots at the top,
-middle, and bottom of the domain (shell).
-"""
-
-def slice(fname:str, field:str,
-          loc: str = "top", width_factor: float = 3.0) -> None:
-    """
-    A slice plot of the dataset for Spherical2D geometry.
-
-    Parameter
-    =======================
-    fname: plot file name
-    field: field parameter
-    loc: location on the domain. {top, mid, bot}
+def slice(fnames:List[str], fields:List[str],
+          loc: str = "top", widthScale: float = 3.0,
+          theta: Optional[float] = None) -> None:
     """
+    A slice plot of the datasets for different field parameters for Spherical2D geometry.
 
-    ds = CastroDataset(fname)
-    currentTime = ds.current_time.in_units("s")
-    print(f"Current time of this plot file is {currentTime} s")
-
-    # Some geometry properties
-    rr = ds.domain_right_edge[0].in_units("km")
-    rl = ds.domain_left_edge[0].in_units("km")
-    dr = rr - rl
-    r_center = 0.5 * (rr + rl)
+    Parameters
+    ==================================================================================
+    fnames: A list of file names to plot multiple slice plots between different
+            plot files for a given field parameter.
+            Note that either fname or field must be single valued.
 
-    thetar = ds.domain_right_edge[1]
-    thetal = ds.domain_left_edge[1]
-    dtheta = thetar - thetal
-    theta_center = 0.5 * (thetar + thetal)
+    fields: A list of field parameters to plot multiple slice plots between different
+            field parameters for a given file.
+            Note that either fname or field must be single valued.
 
-    # Domain width of the slice plot
-    width = width_factor * dr
-    box_widths = (width, width)
+    loc:    preset center location of the domain. {top, mid, bot}
 
-    loc = loc.lower()
-    loc_options = ["top", "mid", "bot"]
-
-    if loc not in loc_options:
-        raise Exception("loc parameter must be top, mid or bot")
+    widthScale: scaling for the domain width of the slice plot
 
-    # Centers for the Top, Mid and Bot panels
-    centers = {"top":(r_center*np.sin(thetal)+0.5*width, r_center*np.cos(thetal)),
-               "mid":(r_center*np.sin(theta_center), r_center*np.cos(theta_center)),
-               "bot":(r_center*np.sin(thetar)+0.5*width, r_center*np.cos(thetar))}
-
-    # Note we can also set center during SlicePlot, however then we would enter in [r_center, theta_center, 0]
-    # rather than the physical R-Z coordinate if we do it via sp.set_center
-
-    sp = yt.SlicePlot(ds, 'phi', field, width=box_widths)
-    sp.set_center(centers[loc])
+    theta:  user defined theta center of the slice plot
+    """
 
-    sp.set_cmap(field, "viridis")
-    if field in ["x_velocity", "y_velocity", "z_velocity"]:
-        sp.set_cmap(field, "coolwarm")
-    elif field == "Temp":
-        sp.set_zlim(field, 5.e7, 2.5e9)
-        sp.set_cmap(field, "magma_r")
-    elif field == "enuc":
-        sp.set_zlim(field, 1.e18, 1.e20)
-    elif field == "density":
-        sp.set_zlim(field, 1.e-3, 5.e8)
+    ts = [CastroDataset(fname) for fname in fnames]
+
+    fig = plt.figure(figsize=(16, 9))
+
+    num = len(fields)*len(fnames)
+    ny = math.ceil(np.sqrt(num))
+    nx = math.ceil(num/ny)
+
+    grid = ImageGrid(fig, 111, nrows_ncols=(nx, ny),
+                     axes_pad=1, label_mode="L", cbar_location="right",
+                     cbar_mode="each", cbar_size="2.5%", cbar_pad="0%")
+
+    # Output plot file name
+    outName = "xrb_spherical_slice.png"
+
+    for j, ds in enumerate(ts):
+        # Process information for each dataset
+
+        # Some geometry properties
+        rr = ds.domain_right_edge[0].in_units("km")
+        rl = ds.domain_left_edge[0].in_units("km")
+        dr = rr - rl
+        r_center = 0.5 * (rr + rl)
+
+        thetar = ds.domain_right_edge[1]
+        thetal = ds.domain_left_edge[1]
+        dtheta = thetar - thetal
+        theta_center = 0.5 * (thetar + thetal)
+
+        # Domain width of the slice plot
+        width = widthScale * dr
+        box_widths = (width, width)
+
+        # Preset centers for the Top, Mid and Bot panels
+        # Centers will be physical coordinates in Cylindrical, i.e. R-Z
+        centers = {"top":(r_center*np.sin(thetal)+0.5*width, r_center*np.cos(thetal)),
+                   "mid":(r_center*np.sin(theta_center), r_center*np.cos(theta_center)),
+                   "bot":(r_center*np.sin(thetar)+0.5*width, r_center*np.cos(thetar))}
+
+        if theta is None:
+            center = centers[loc]
+        else:
+            R = r_center*np.sin(theta)
+            Z = r_center*np.cos(theta)
+            if R < 0.5*width:
+                R = 0.5*width
+            center = [R, Z]
+
+        for i, field in enumerate(fields):
+            # Plot each field parameter
+
+            sp = yt.SlicePlot(ds, 'phi', field, width=box_widths, fontsize=16)
+            sp.set_center(center)
+
+            sp.set_cmap(field, "viridis")
+            if field in ["x_velocity", "y_velocity", "z_velocity"]:
+                sp.set_cmap(field, "coolwarm")
+                if field == "z_velocity":
+                    sp.set_zlim(field, -2.e8, 2.e8)
+                    sp.set_log(field, False)
+            elif field == "Temp":
+                sp.set_zlim(field, 5.e7, 2.5e9)
+                sp.set_cmap(field, "magma_r")
+            elif field == "enuc":
+                sp.set_zlim(field, 1.e18, 1.e20)
+            elif field == "density":
+                sp.set_zlim(field, 1.e-3, 5.e8)
+
+            sp.set_buff_size((2400,2400))
+            sp.set_axes_unit("km")
+            # sp.annotate_text((0.05, 0.05), f"{currentTime.in_cgs():8.5f} s")
+
+            plot = sp.plots[field]
+            plot.figure = fig
+            plot.axes = grid[i+j*len(fields)].axes
+            plot.cax = grid.cbar_axes[i+j*len(fields)]
+
+            sp._setup_plots()
+
+    if len(fnames) == 1:
+        time = ts[0].current_time.in_units("ms")
+        if float(time) < 1e-1:
+            time = ts[0].current_time.in_units("us")
+
+        # Determine position of the text on grid
+        xyPositions = {(1, 1): (0.78, 0.02),
+                       (1, 2): (0.95, 0.075),
+                       (2, 2): (0.78, 0.02),
+                       (2, 3): (0.9, 0.02),
+                       (3, 3): (0.78, 0.02)
+                       }
+        xPosition, yPosition = xyPositions.get((nx, ny), (0.78, 0.02))
+
+        fig.text(xPosition, yPosition, f"t = {time:.3f}", fontsize=16,
+                 horizontalalignment='right', verticalalignment='bottom',
+                 color="black", transform=fig.transFigure)
+
+        outName = f"{ts[0]}_slice.png"
+
+    fig.set_size_inches(16, 9)
+    fig.tight_layout()
+    fig.savefig(outName, format="png", bbox_inches="tight")
 
-    # sp.annotate_text((0.05, 0.05), f"{currentTime.in_cgs():8.5f} s")
-    sp.save(f"{ds}_{loc}")
 
 if __name__ == "__main__":
 
-    if len(sys.argv) < 3:
-        raise Exception("Please enter parameters in order of: fname field_name width_factor[optional] loc[optional]")
-
-    fname = sys.argv[1]
-    field = sys.argv[2]
-    loc = "top"
-    width_factor = 3.0
+    parser = argparse.ArgumentParser(description="""
+        A slice plot script for xrb_spherical problem.
+        Given a list of plotfiles or a list of field parameters,
+        it plots multiple slice plots.
+        """)
+
+    parser.add_argument('--fnames', nargs='+', type=str,
+                        help="""dataset file names for plotting. Accepts one or more datasets.
+                        If multiple file names are given, a grid of slice plots of different
+                        files will be plotted for a given field parameter.
+                        Note that either fnames or field must be single valued.""")
+    parser.add_argument('--fields', nargs='+', type=str,
+                        help="""field parameters for plotting. Accepts one or more datasets.
+                        If multiple parameters are given, a grid of slice plots of different
+                        field parameters will be plotted for a given fname.
+                        Note that either fnames or fields must be single valued.
+                        """)
+    parser.add_argument('-l', '--loc', default='top', type=str, metavar="{top, mid, bot}",
+                        help="""preset center location of the plot domain.
+                        Enter one of the three choices: {top, mid, bot}""")
+    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('-w', '--width', default=4.0, type=float,
+                        help="scaling for the domain width of the slice plot")
+
+
+    args = parser.parse_args()
+
+    if len(args.fnames) > 1 and len(args.fields) > 1:
+        parser.error("Either fnames or fields must be single valued!")
+
+    loc = args.loc.lower()
+    loc_options = ["top", "mid", "bot"]
 
-    if len(sys.argv) == 4:
-        width_factor = float(sys.argv[3])
-    elif len(sys.argv) > 4:
-        loc = sys.argv[4]
+    if loc not in loc_options:
+        parser.error("loc must be one of the three: {top, mid, bot}")
 
-    slice(fname, field, loc=loc, width_factor=width_factor)
+    slice(args.fnames, args.fields,
+          loc=loc, theta=args.theta, widthScale=args.width)

Exec/science/xrb_spherical/inputs.He.1000Hz

@@ -6,7 +6,7 @@ stop_time =  3.0
 geometry.is_periodic = 0       0
 geometry.coord_sys   = 2                  # 0 => cart, 1 => RZ  2=>spherical
 geometry.prob_lo     = 1.1e6           0
-geometry.prob_hi     = 1.13072e6       3.1415926
+geometry.prob_hi     = 1.13072e6       3.141592653589793
 amr.n_cell           = 768            2304  #192             1152
 
 # >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
@@ -95,16 +95,16 @@ amr.max_grid_size   = 128
 amr.n_error_buf     = 2 2 2 2 # number of buffer cells in error est
 
 # CHECKPOINT FILES
-amr.check_file      = flame_wave_1000Hz_chk        # root name of checkpoint file
+amr.check_file      = xrb_spherical_1000Hz_chk        # root name of checkpoint file
 amr.check_int       = 250        # number of timesteps between checkpoints
 
 # PLOTFILES
-amr.plot_file        = flame_wave_1000Hz_plt        # root name of plotfile
+amr.plot_file        = xrb_spherical_1000Hz_plt        # root name of plotfile
 amr.plot_per         = 5.e-3      # number of seconds between plotfiles
 amr.derive_plot_vars = ALL
 
-amr.small_plot_file        = flame_wave_1000Hz_smallplt        # root name of plotfile
-amr.small_plot_per         = 1.e-6 #1.e-4      # number of seconds between plotfiles
+amr.small_plot_file         = xrb_spherical_1000Hz_smallplt        # root name of plotfile
+amr.small_plot_per          = 1.e-6 #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