From 26def2f432854a89bf84c2210df930c94e1b7d45 Mon Sep 17 00:00:00 2001 From: zhi Date: Sun, 22 Dec 2024 17:36:18 -0500 Subject: [PATCH] update script, makefile and input --- Exec/science/xrb_spherical/GNUmakefile | 4 +- Exec/science/xrb_spherical/analysis/slice.py | 235 +++++++++++++------ Exec/science/xrb_spherical/inputs.He.1000Hz | 10 +- 3 files changed, 173 insertions(+), 76 deletions(-) diff --git a/Exec/science/xrb_spherical/GNUmakefile b/Exec/science/xrb_spherical/GNUmakefile index d53475be8d..f77c92cf66 100644 --- a/Exec/science/xrb_spherical/GNUmakefile +++ b/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 diff --git a/Exec/science/xrb_spherical/analysis/slice.py b/Exec/science/xrb_spherical/analysis/slice.py index c6a9649160..40d202e16b 100755 --- a/Exec/science/xrb_spherical/analysis/slice.py +++ b/Exec/science/xrb_spherical/analysis/slice.py @@ -3,92 +3,189 @@ import sys import os import yt +import argparse +from typing import List, Optional import numpy as np import matplotlib.pyplot as plt 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: +SMALL_SIZE = 28 +MEDIUM_SIZE = 30 +BIGGER_SIZE = 32 + +plt.rc('font', size=SMALL_SIZE) # controls default text sizes +plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title +plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels +plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels +plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels +plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize +plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title +plt.rc('xtick.major', size=7, width=2) +plt.rc('xtick.minor', size=5, width=1) +plt.rc('ytick.major', size=7, width=2) +plt.rc('ytick.minor', size=5, width=1) + +def slice(fnames:List[str], fields:List[str], + loc: str = "top", widthScale: float = 4.0, + coord: Optional[List[float, float]] = None) -> 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} - """ - - 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]) + coord: user defined center of the slice in the format [r_center, theta_center] + """ - 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(fnames) + fig = plt.figures(figsize=(16, 9)) + grid = ImageGrid(fig, 111, nrows_ncols=(len(fields)*len(ts), 1), + axes_pad=0.25, label_mode="L", cbar_location="right", + cbar_mode="each", cbar_size="2.5%", cbar_pad="0%") + + # We will check if all files have the same timestamp later. + # This matters on how we annotate timestamp for the slice plots. + if len(fnames) == 1: + sameTime = True + else: + refTimeStamp = ts[0].current_time.in_units("ms") + sameTime = all(refTimeStamp == ds.current_time.in_units("ms") for ds in ts) + + for j, ds in enumerate(ts): + # Process information for each dataset + + currentTime = ds.current_time.in_units("ms") + print(f"Current time of this plot file is {currentTime}") + + # 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 = {"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 coord is None: + # Set center vis sp.set_center. + #This will be physical coordinates in Cartesian + center = centers[loc] + else: + # Set center during SlicePlot call. + # This will be in Spherical coordinate: [r_center, theta_center, 0] + # coord will be in [r_center, theta_center], so append 0 here. + center = coord.append(0) + + for i, field in enumerate(fields): + # Plot each field parameter + + # 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 + + if coord is None: + sp = yt.SlicePlot(ds, 'phi', field, width=box_widths) + sp.set_center(center) + else: + sp = yt.SlicePlot(ds, 'phi', field, center=center, width=box_widths) + + 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) + + 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)] + + if i+j*len(fields) < len(fnames)+len(fields)-1: + grid[i].axes.xaxis.offsetText.set_visible(False) + + sp._setup_plots() + + if sameTime: + fig.text(0.8, 0.05, f"t = {ds.current_time}", + horizontalalignment='right', verticalalignment='center', + color="black", transform=fig.transFigure) + + fig.set_size_inches(24, 13.5) + fig.tight_layout() + fig.savefig(f"{ds}",format="pdf",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 plot file and field name, it gives slice plots at the top, + middle, and bottom of the domain (shell). + """) + + 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, + help="""preset center location of the plot domain. + Enter one of the three choices: {top, mid, bot}""") + parser.add_argument('-c', '--coord', nargs=2, type=float, + help="""user defined center location of the plot domain. + Enter two numbers in the format of [r_center, theta_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, coord=args.coord, width_factor=args.width) diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz index 17864ed3f3..b2fa2b5fe0 100644 --- a/Exec/science/xrb_spherical/inputs.He.1000Hz +++ b/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