sync up with the MAESTROeX xrb_layered (#2830)

also add a non-uniform grid version

Exec/science/xrb_layered/_prob_params

@@ -1,4 +1,8 @@
-H_min                 real           1.e-4_rt      y
+# H mass fraction above which to refine
+H_min                 real           -1.0          y
+
+# He mass fraction above which to refine
+He_min                real           -1.0          y
 
 cutoff_density        real           50.0_rt       y
 

Exec/science/xrb_layered/analysis/slice_vertical.py

@@ -0,0 +1,113 @@
+import argparse
+import os
+import sys
+
+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 amu, cm
+
+matplotlib.use('agg')
+
+def doit(plotfiles):
+
+    plotfile = sys.argv[1]
+    ds = CastroDataset(plotfile)
+
+    xmin = ds.domain_left_edge[0]
+    xmax = ds.domain_right_edge[0]
+    xctr = 0.5 * (xmin + xmax)
+    L_x = xmax - xmin
+
+    ymin = ds.domain_left_edge[1]
+    ymax = ds.domain_right_edge[1]
+    yctr = 0.5*(ymin + ymax)
+
+    zmin = ds.domain_left_edge[2]
+    zmax = ds.domain_right_edge[2]
+    zctr = 0.5*(zmin + zmax)
+
+
+    zctr = 0.5*(zmin + zmax)
+    L_z = zmax - zmin
+
+    fig = plt.figure()
+
+    grid = ImageGrid(fig, 111, nrows_ncols=(len(plotfiles), 1),
+                     axes_pad=0.25, label_mode="L",
+                     cbar_mode="single", cbar_size="3%")
+
+    for i, pf in enumerate(plotfiles):
+
+        ds = CastroDataset(pf)
+
+        f = "magvort"
+
+        sp = yt.SlicePlot(ds, "y", f, origin="native", center=[xctr, yctr, zctr],
+                          width=[L_z, L_x], fontsize="11")
+        sp.set_buff_size((4800,4800))
+        sp.swap_axes()
+
+        if f == "Temp":
+            sp.set_zlim(f, 5.e5, 5.e8)
+            sp.set_cmap(f, "magma_r")
+        elif f == "enuc":
+            sp.set_zlim(f, 1.e14, 1.e17)
+            sp.set_log(f, True)
+            sp.set_cmap(f, "plasma_r")
+        elif f == "density":
+            sp.set_zlim(f, 1.e-5, 5.e6)
+        elif f == "abar":
+            sp.set_zlim(f, 1, 8)
+            sp.set_log(f, False)
+            sp.set_cmap(f, "plasma_r")
+        elif f == "magvort":
+            sp.set_zlim(f, 2.0, 2.e5)
+            sp.set_log(f, True)
+        elif f.startswith(r"X("):
+            sp.set_zlim(f, 1.e-5, 1.0)
+            sp.set_log(f, True)
+
+        sp.set_axes_unit("cm")
+
+        sp.annotate_text((0.05, 0.05), "{:5.2f} ms".format(1000.0*float(ds.current_time.in_cgs())),
+                         coord_system="figure", text_args={"color": "black", "size": 11})
+
+        plot = sp.plots[f]
+        plot.figure = fig
+        plot.axes = grid[i].axes
+        plot.cax = grid.cbar_axes[i]
+        if i < len(plotfiles)-1:
+            grid[i].axes.xaxis.offsetText.set_visible(False)
+
+        sp._setup_plots()
+
+    fig.set_size_inches(8, 6.5)
+    fig.subplots_adjust(left=0.05, right=0.95, top=0.95, bottom=0.05)
+    fig.tight_layout()
+    fig.savefig(f"xrb_{f}.png")
+
+if __name__ == "__main__":
+
+    p = argparse.ArgumentParser()
+
+    p.add_argument("--skip", type=int, default=1,
+                   help="interval between plotfiles")
+    p.add_argument("plotfiles", type=str, nargs="+",
+                   help="list of plotfiles to plot")
+
+    args = p.parse_args()
+
+    plot_prefix = args.plotfiles[0].split("plt")[0] + "plt"
+    plot_nums = sorted([p.split("plt")[1] for p in args.plotfiles], key=int)
+
+    plotfiles = []
+    for n in range(0, len(plot_nums), args.skip):
+        plotfiles.append("{}{}".format(plot_prefix, plot_nums[n]))
+
+    doit(plotfiles)

Exec/science/xrb_layered/analysis/vol-zvel.py

@@ -0,0 +1,99 @@
+#!/usr/bin/env python
+
+import sys
+
+import matplotlib
+import numpy as np
+
+import yt
+from yt.frontends.boxlib.api import CastroDataset
+from yt.units import cm
+#from yt.visualization.volume_rendering.render_source import VolumeSource
+from yt.visualization.volume_rendering.api import Scene, create_volume_source
+
+matplotlib.use('agg')
+
+resolution = (1920, 1080)
+
+# this is for the wdconvect problem
+
+def doit(plotfile):
+
+    ds = CastroDataset(plotfile)
+    ds._periodicity = (True, True, True)
+
+    field = ('boxlib', 'z_velocity')
+    ds._get_field_info(field).take_log = False
+
+    sc = Scene()
+
+
+    # add a volume: select a sphere
+    #center = (0, 0, 0)
+    #R = (5.e8, 'cm')
+
+    #dd = ds.sphere(center, R)
+
+    vol = create_volume_source(ds.all_data(), field=field)
+    sc.add_source(vol)
+
+
+    # transfer function
+    vals = [-1.e7, -5.e6, -2.5e6, 2.5e6, 5.e6, 1.e7]
+    sigma = 5.e5
+
+    tf =  yt.ColorTransferFunction((min(vals), max(vals)))
+
+    tf.clear()
+
+    cmap = "coolwarm"
+
+    for v in vals:
+        tf.sample_colormap(v, sigma**2, alpha=0.2, colormap=cmap)
+
+    sc.get_source(0).transfer_function = tf
+
+    cam = sc.add_camera(ds, lens_type="perspective")
+    cam.resolution = resolution
+
+    # view 1
+
+    center = 0.5*(ds.domain_left_edge + ds.domain_right_edge)
+
+    cam.position = [2.5*ds.domain_right_edge[0],
+                    2.5*ds.domain_right_edge[1],
+                    center[2]+0.25*ds.domain_right_edge[2]]
+
+    # look toward the center -- we are dealing with an octant
+    normal = (center - cam.position)
+    normal /= np.sqrt(normal.dot(normal))
+
+    cam.switch_orientation(normal_vector=normal,
+                           north_vector=[0., 0., 1.])
+    cam.set_width(ds.domain_width)
+
+    sc.camera = cam
+
+    sc.save(f"{plotfile}_zvel.png", sigma_clip=4.0)
+
+    sc.annotate_axes(alpha=0.005, thickness=6)
+    sc.annotate_domain(ds, color=np.array([0.05, 0.05, 0.05, 0.05]))
+    sc.save(f"{plotfile}_zvel_axes.png", sigma_clip=3.0)
+
+    sc.save_annotated(f"{plotfile}_zvel_annotated.png",
+                      sigma_clip=4.0, label_fmt="%.2f",
+                      text_annotate=[[(0.05, 0.05),
+                                      f"t = {ds.current_time.d:7.5f} s",
+                                      dict(horizontalalignment="left")]])
+
+
+if __name__ == "__main__":
+
+    # Choose a field
+    plotfile = ""
+
+
+    try: plotfile = sys.argv[1]
+    except: sys.exit("ERROR: no plotfile specified")
+
+    doit(plotfile)

Exec/science/xrb_layered/inputs_2d

@@ -7,7 +7,7 @@ geometry.is_periodic = 1       0
 geometry.coord_sys   = 0                  # 0 => cart, 1 => RZ  2=>spherical
 geometry.prob_lo     = 0.0     0.0
 geometry.prob_hi     = 3072    3072
-amr.n_cell           = 128     128
+amr.n_cell           = 256     256
 
 # >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
 # 0 = Interior           3 = Symmetry
@@ -75,11 +75,14 @@ amr.plot_int        = 1000       # number of timesteps between plotfiles
 amr.derive_plot_vars = ALL
 
 # PROBLEM PARAMETERS
-problem.model_name =  "toy_atm_hot_castro_3cm.hse"
+problem.model_name =  "toy_atm_hot_3cm.hse"
 
 problem.apply_perturbation = 1
 
-problem.model_shift = 300.0
+problem.H_min = 1.e-5
+problem.He_min = 1.e-5
+
+problem.model_shift = 750.0
 
 # MICROPHYSICS
 integrator.jacobian = 1

Exec/science/xrb_layered/inputs_2d.nonsquare

@@ -0,0 +1,96 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 1000000
+stop_time =  1000.0
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic = 1       0
+geometry.coord_sys   = 0                  # 0 => cart, 1 => RZ  2=>spherical
+geometry.prob_lo     = 0.0     0.0
+geometry.prob_hi     = 3072    3072
+amr.n_cell           = 128     256
+
+castro.allow_non_unit_aspect_zones = 1
+castro.domain_is_plane_parallel = 1
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+castro.lo_bc       =  0   3
+castro.hi_bc       =  0   2
+
+castro.fill_ambient_bc = 1
+castro.ambient_fill_dir = 1
+castro.ambient_outflow_vel = 1
+
+
+# WHICH PHYSICS
+castro.do_hydro = 1
+castro.do_react = 1
+castro.add_ext_src   = 0
+castro.do_grav = 1
+castro.do_sponge = 1
+
+castro.ppm_type = 1
+castro.grav_source_type = 2
+castro.use_pslope = 1
+castro.pslope_cutoff_density = 1.e4
+
+gravity.gravity_type = ConstantGrav
+gravity.const_grav   = -1.29e14
+
+# burning
+castro.react_rho_min = 1.e2
+castro.react_T_min = 5.e6
+
+
+# TIME STEP CONTROL
+castro.cfl            = 0.8     # cfl number for hyperbolic system
+castro.init_shrink    = 0.1     # scale back initial timestep
+castro.change_max     = 1.1     # max time step growth
+
+# SPONGE
+castro.sponge_upper_density = 1.e3
+castro.sponge_lower_density = 1.e2
+castro.sponge_timescale     = 1.0e-7
+
+# DIAGNOSTICS & VERBOSITY
+castro.sum_interval   = 1       # timesteps between computing mass
+castro.v              = 1       # verbosity in Castro.cpp
+amr.v                 = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING 
+amr.max_level       = 1       # maximum level number allowed
+amr.ref_ratio       = 2 2 2 2 # refinement ratio
+amr.regrid_int      = 2 2 2 2 # how often to regrid
+amr.blocking_factor = 8       # block factor in grid generation
+amr.max_grid_size   = 256
+amr.n_error_buf     = 2 2 2 2 # number of buffer cells in error est
+
+# CHECKPOINT FILES
+amr.check_file      = xrb_chk    # root name of checkpoint file
+amr.check_int       = 1000       # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_file       = xrb_plt    # root name of plotfile
+amr.plot_per        = 0.001
+amr.derive_plot_vars = ALL
+
+# PROBLEM PARAMETERS
+problem.model_name =  "toy_atm_hot_3cm.hse"
+
+problem.apply_perturbation = 1
+
+problem.model_shift = 750.0
+
+problem.H_min = -1
+problem.He_min = -1
+
+# MICROPHYSICS
+integrator.jacobian = 1
+
+integrator.atol_spec = 1.e-6
+integrator.rtol_spec = 1.e-6
+
+integrator.X_reject_buffer = 4.0