xrb spherical problem setup (#2972)

Initial problem setup for full-star xrb flame in spherical shell based on the flame_wave problem.

Exec/science/xrb_spherical/GNUmakefile

@@ -0,0 +1,40 @@
+PRECISION  = DOUBLE
+PROFILE    = FALSE
+
+DEBUG      = FALSE
+
+DIM        = 2
+
+COMP       = gnu
+
+USE_MPI    = TRUE
+
+USE_GRAV   = TRUE
+USE_REACT  = FALSE
+
+USE_ROTATION  = FALSE
+USE_DIFFUSION = FALSE
+
+# define the location of the CASTRO top directory
+CASTRO_HOME  ?= ../../..
+
+USE_JACOBIAN_CACHING = TRUE
+USE_MODEL_PARSER = TRUE
+NUM_MODELS := 2
+
+# This sets the EOS directory in $(MICROPHYSICS_HOME)/eos
+EOS_DIR     := helmholtz
+
+# This sets the network directory in $(MICROPHYSICS_HOME)/networks
+NETWORK_DIR := subch_base
+
+INTEGRATOR_DIR := VODE
+
+CONDUCTIVITY_DIR := stellar
+
+PROBLEM_DIR ?= ./
+
+Bpack   := $(PROBLEM_DIR)/Make.package
+Blocs   := $(PROBLEM_DIR)
+
+include $(CASTRO_HOME)/Exec/Make.Castro

Exec/science/xrb_spherical/Make.package

@@ -0,0 +1,2 @@
+CEXE_headers += initial_model.H
+

Exec/science/xrb_spherical/README.md

@@ -0,0 +1,6 @@
+# xrb_spherical
+
+This is the full-star XRB flame setup based on flame_wave.
+This setup uses a spherical 2D geometry to model XRB flame
+on a spherical shell with initial temperature perturbation
+on the north pole.

Exec/science/xrb_spherical/_prob_params

@@ -0,0 +1,68 @@
+
+dtemp                   real              3.81e8_rt                 y
+
+theta_half_max          real              1.745e-2_rt               y
+
+theta_half_width        real              4.9e-3_rt                 y
+
+# cutoff mass fraction of the first species for refinement
+X_min                   real              1.e-4_rt                  y
+
+# do we dynamically refine based on density? or based on height?
+tag_by_density          integer           1                         y
+
+# used for tagging if tag_by_density = 1
+cutoff_density          real              500.e0_rt                 y
+
+# used if we are refining based on height rather than density
+refine_height           real              3600                      y
+
+T_hi                    real              5.e8_rt                   y
+
+T_star                  real              1.e8_rt                   y
+
+T_lo                    real              5.e7_rt                   y
+
+dens_base               real              2.e6_rt                   y
+
+H_star                  real              500.e0_rt                 y
+
+atm_delta               real              25.e0_rt                  y
+
+fuel1_name              string         "helium-4"                y
+
+fuel2_name              string         ""                        y
+
+fuel3_name              string         ""                        y
+
+fuel4_name              string         ""                        y
+
+ash1_name               string         "iron-56"                 y
+
+ash2_name               string         ""                        y
+
+ash3_name               string         ""                        y
+
+fuel1_frac              real              1.0_rt                    y
+
+fuel2_frac              real              0.0_rt                    y
+
+fuel3_frac              real              0.0_rt                    y
+
+fuel4_frac              real              0.0_rt                    y
+
+ash1_frac               real              1.0_rt                    y
+
+ash2_frac               real              0.0_rt                    y
+
+ash3_frac               real              0.0_rt                    y
+
+low_density_cutoff      real              1.e-4_rt                  y
+
+smallx                  real              1.e-10_rt                 y
+
+r_refine_distance       real              1.e30_rt                  y
+
+max_hse_tagging_level   integer           2                         y
+
+max_base_tagging_level  integer           1                         y

Exec/science/xrb_spherical/analysis/r_profile.py

@@ -0,0 +1,56 @@
+#!/usr/bin/env python3
+
+# Spherical R profile at different theta
+
+import os
+import sys
+import yt
+import matplotlib.pyplot as plt
+import numpy as np
+from functools import reduce
+import itertools
+
+import matplotlib.ticker as ptick
+from yt.frontends.boxlib.api import CastroDataset
+from yt.units import cm
+
+
+plotfile = sys.argv[1]
+ds = CastroDataset(plotfile)
+
+rmin = ds.domain_left_edge[0]
+rmax = rmin + 5000.0*cm
+#rmax = ds.domain_right_edge[0]
+print(ds.domain_left_edge[1])
+fig, _ax = plt.subplots(2,2)
+
+axes = list(itertools.chain(*_ax))
+
+fig.set_size_inches(7.0, 8.0)
+
+fields = ["Temp", "density", "x_velocity", "y_velocity"]
+nice_names = [r"$T$ (K)", r"$\rho$ (g/${cm}^3$)", r"$u$ (cm/s)", r"$v$ (cm/s)"]
+
+# 4 rays at different theta values
+thetal = ds.domain_left_edge[1]
+thetar = ds.domain_right_edge[1]
+thetas = [thetal, 0.25*thetar, 0.5*thetar, 0.75*thetar]
+
+for i, f in enumerate(fields):
+
+    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))
+        isrt = np.argsort(ray["t"])
+        axes[i].plot(ray['r'][isrt], ray[f][isrt], label=r"$\theta$ = {:.4f}".format(float(theta)))
+
+    axes[i].set_xlabel(r"$r$ (cm)")
+    axes[i].set_ylabel(nice_names[i])
+    axes[i].set_yscale("symlog")
+
+    if i == 0:
+        axes[0].legend(frameon=False, loc="lower left")
+
+#fig.set_size_inches(10.0, 9.0)
+plt.tight_layout()
+plt.savefig("{}_profiles.png".format(os.path.basename(plotfile)))

Exec/science/xrb_spherical/analysis/slice.py

@@ -0,0 +1,94 @@
+#!/usr/bin/env python3
+
+import sys
+import os
+import yt
+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:
+    """
+    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)
+
+    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 = width_factor * dr
+    box_widths = (width, width)
+
+    loc = loc.lower()
+    loc_options = ["top", "mid", "bot"]
+
+    if loc not in loc_options:
+        raise Exception("loc parameter must be top, mid or bot")
+
+    # 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])
+
+    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(f, 5.e7, 2.5e9)
+        sp.set_cmap(f, "magma_r")
+    elif field == "enuc":
+        sp.set_zlim(f, 1.e18, 1.e20)
+    elif field == "density":
+        sp.set_zlim(f, 1.e-3, 5.e8)
+
+    # 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
+
+    if len(sys.argv) == 4:
+        width_factor = float(sys.argv[3])
+    elif len(sys.argv) > 4:
+        loc = sys.argv[4]
+
+    slice(fname, field, loc=loc, width_factor=width_factor)

Exec/science/xrb_spherical/initial_model.H

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