Skip to content

Commit

Permalink
update script, makefile and input
Browse files Browse the repository at this point in the history
  • Loading branch information
@zhichen3
zhichen3 committed Dec 22, 2024
1 parent 34d2b9a commit 26def2f
Show file tree
Hide file tree
Showing 3 changed files with 173 additions and 76 deletions.
4 changes: 2 additions & 2 deletions Exec/science/xrb_spherical/GNUmakefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
235 changes: 166 additions & 69 deletions Exec/science/xrb_spherical/analysis/slice.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
10 changes: 5 additions & 5 deletions Exec/science/xrb_spherical/inputs.He.1000Hz
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <<<<<<<<<<<<<<<<
Expand Down Expand Up @@ -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

Expand Down

0 comments on commit 26def2f

Please sign in to comment.