a setup for exploring shock burning

Exec/science/Detonation/shock_paper/README.md

@@ -0,0 +1,23 @@
+# Shock Burning Experiments
+
+This directory is meant to explore shock burning with detonation.  Compile as:
+
+```
+make USE_SIMPLIFIED_SDC=TRUE USE_SHOCK_VAR=TRUE NETWORK_DIR=aprox13 -j 4
+```
+
+Then setup a suite of simulations with the following resolutions (the
+`inputs-det-x.subch_base`) here is for the coarsest run:
+
+
+| resolution   |  base grid  |  levels (4x jumps)  |
+| ------------ | ----------- | ------------------- |
+|   12.288 km  |       48    |         1           |
+|    1.536 km  |      384    |         1           |
+|    0.192 km  |     3072    |         1           |
+|     2400 cm  |     6144    |         2           |
+|      300 cm  |    12288    |         3           |
+
+
+
+

Exec/science/Detonation/shock_paper/det_speed_comp.py

@@ -0,0 +1,31 @@
+import matplotlib.pyplot as plt
+
+import detonation
+
+runs = [("res12.288km", 12.288),
+        ("res1.536km", 1.536),
+        ("res0.192km", 0.192),
+        ("res0.024km", 0.024),
+        ("res0.003km", 0.003)]
+
+res = []
+v = []
+dv = []
+
+for ddir, dx in runs:
+    res.append(dx)
+    d = detonation.Detonation(ddir)
+    v.append(d.v)
+    dv.append(d.v_sigma)
+
+
+fig, ax = plt.subplots()
+
+ax.errorbar(res, v, yerr=dv, fmt="o")
+
+ax.set_xlabel(r"$\Delta x$ (km)")
+ax.set_ylabel(r"$v_\mathrm{det}$ (cm / s)")
+
+ax.set_xscale("log")
+
+fig.savefig("det_speeds.png")

Exec/science/Detonation/shock_paper/detonation.py

@@ -0,0 +1,105 @@
+import glob
+import os
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import yt
+from yt.frontends.boxlib.api import CastroDataset
+
+
+yt.funcs.mylog.setLevel(50)
+
+class Profile:
+    """read a plotfile using yt and store the 1d profile for T and enuc"""
+
+    def __init__(self, plotfile):
+        ds = CastroDataset(plotfile)
+
+        time = float(ds.current_time)
+        ad = ds.all_data()
+
+        # Sort the ray values by 'x' so there are no discontinuities
+        # in the line plot
+        srt = np.argsort(ad['x'])
+        x_coord = np.array(ad['x'][srt])
+        temp = np.array(ad['Temp'][srt])
+        enuc = np.array(ad['enuc'][srt])
+        shock = np.array(ad['Shock'][srt])
+
+        self.time = time
+        self.x = x_coord
+        self.T = temp
+        self.enuc = enuc
+        self.shock = shock
+
+    def find_x_for_T(self, T_0=1.e9):
+        """ given a profile x(T), find the x_0 that corresponds to T_0 """
+
+        # our strategy here assumes that the hot ash is in the early
+        # part of the profile.  We then find the index of the first
+        # point where T drops below T_0
+        try:
+            idx = np.where(self.T < T_0)[0][0]
+        except IndexError:
+            idx = len(self.T)-1
+
+        T1 = self.T[idx-1]
+        x1 = self.x[idx-1]
+
+        T2 = self.T[idx]
+        x2 = self.x[idx]
+
+        slope = (x2 - x1)/(T2 - T1)
+
+        return x1 + slope*(T_0 - T1)
+
+
+class Detonation:
+    def __init__(self, dirname):
+        self.name = dirname
+
+        self.v = None
+        self.v_sigma = None
+
+        # find all the output (plot) files
+        self.files = glob.glob(f"{dirname}/*plt?????")
+        self.files.sort()
+
+        # precompute the velocity and the data profiles
+        if len(self.files) >= 3:
+            self.v, self.v_sigma = self.get_velocity()
+        else:
+            self.v, self.v_sigma = 0.0, 0.0
+
+    def get_velocity(self):
+        """look at the last 2 plotfiles and estimate the velocity by
+        finite-differencing"""
+
+        vs = []
+        pairs = [(-2, -1), (-3, -1), (-3, -2)]
+
+        for i1, i2 in pairs:
+            p1 = self.get_data(i1)
+            p2 = self.get_data(i2)
+
+            # we'll do this by looking at 3 different temperature
+            # thresholds and averaging
+            T_ref = [1.e9, 2.e9, 3.e9]
+
+            for T0 in T_ref:
+                x1 = p1.find_x_for_T(T0)
+                x2 = p2.find_x_for_T(T0)
+                vs.append((x1 - x2)/(p1.time - p2.time))
+
+        vs = np.array(vs)
+        v = np.mean(vs)
+        v_sigma = np.std(vs)
+        return v, v_sigma
+
+    def get_data(self, num=-1):
+        """get the temperature and energy generation rate from the
+        num'th plotfile (defaults to the last)"""
+
+        return Profile(os.path.join(self.files[num]))
+

Exec/science/Detonation/shock_paper/inputs-det-x.subch_base

@@ -0,0 +1,112 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 200000
+stop_time =  0.02
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic = 0 0 0
+geometry.coord_sys   = 0  # 0 => cart, 1 => RZ  2=>spherical
+geometry.prob_lo     = 0     0     0
+geometry.prob_hi     = 5.89824e7  2500  2500
+amr.n_cell           = 48    16    16
+
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+castro.lo_bc       =  3   4   4
+castro.hi_bc       =  2   4   4
+
+# WHICH PHYSICS
+castro.do_hydro = 1
+castro.do_react = 1
+
+castro.react_rho_min = 100.0
+castro.react_T_min = 5.e7
+
+castro.ppm_type = 1
+castro.ppm_temp_fix = 0
+
+castro.time_integration_method = 3
+
+# castro.transverse_reset_density = 1
+
+castro.small_dens   = 1.e-5
+castro.small_temp     = 1.e7
+
+castro.use_flattening = 1
+
+castro.riemann_solver = 1
+
+# TIME STEP CONTROL
+castro.cfl            = 0.5     # cfl number for hyperbolic system
+castro.init_shrink    = 0.1     # scale back initial timestep
+castro.change_max     = 1.05    # scale back initial timestep
+
+#castro.dtnuc_e = 0.1
+
+# DIAGNOSTICS & VERBOSITY
+castro.sum_interval   = 1       # timesteps between computing mass
+castro.v              = 1       # verbosity in Castro.cpp
+amr.v                 = 1       # verbosity in Amr.cpp
+#amr.grid_log        = grdlog  # name of grid logging file
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0       # 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   = 512
+amr.n_error_buf     = 8 8 8 4 4 2 2 2 2 # number of buffer cells in error est
+
+# CHECKPOINT FILES
+amr.check_file      = det_x_chk  # root name of checkpoint file
+amr.check_int       = 10000         # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_file       = det_x_plt  # root name of plotfile
+amr.plot_per = 2.e-3
+amr.derive_plot_vars = ALL
+
+# problem initialization
+
+problem.T_l = 1.1e9
+problem.T_r = 1.75e8
+
+problem.dens = 3.0e6
+problem.cfrac = 0.0
+problem.nfrac = 0.01
+
+problem.smallx = 1.e-10
+
+problem.idir = 1
+
+problem.w_T = 1.e-3
+problem.center_T = 0.2
+
+# refinement
+
+amr.refinement_indicators = temperr tempgrad
+
+amr.refine.temperr.max_level = 0
+amr.refine.temperr.value_greater = 4.e9
+amr.refine.temperr.field_name = Temp
+
+amr.refine.tempgrad.max_level = 5
+amr.refine.tempgrad.gradient = 1.e8
+amr.refine.tempgrad.field_name = Temp
+
+# Microphysics
+
+network.small_x = 1.e-10
+integrator.SMALL_X_SAFE = 1.e-10
+
+integrator.rtol_spec = 1.e-5
+integrator.atol_spec = 1.e-5
+integrator.rtol_enuc = 1.e-5
+integrator.atol_enuc = 1.e-5
+integrator.jacobian = 1
+
+integrator.X_reject_buffer = 4.0
+

Exec/science/Detonation/shock_paper/profile_compare.py

@@ -0,0 +1,60 @@
+#!/usr/bin/env python3
+
+# Take a sequence of plotfiles and plot T and enuc vs. position
+
+import matplotlib
+import numpy as np
+
+matplotlib.use('agg')
+
+import matplotlib.pyplot as plt
+
+import matplotlib.ticker as mticker
+
+import detonation
+
+def plot_Te(data):
+
+    f = plt.figure()
+
+    f.set_size_inches(7.5, 9.0)
+
+    ax_T = f.add_subplot(211)
+    ax_e = f.add_subplot(212)
+
+    for n, _d in enumerate(data):
+
+        ddir, label = _d
+
+        d = detonation.Detonation(ddir)
+        profile = d.get_data(-1)
+
+        ax_T.plot(profile.x, profile.T, label=label)
+        ax_e.plot(profile.x, profile.enuc, label=label)
+
+    ax_T.set_ylabel("T (K)")
+    ax_T.yaxis.set_major_formatter(mticker.ScalarFormatter(useMathText=True))
+    ax_T.xaxis.set_major_formatter(mticker.ScalarFormatter(useMathText=True))
+    ax_T.legend()
+
+    ax_e.set_yscale("log")
+    ax_e.set_ylabel(r"$S_\mathrm{nuc}$ (erg/g/s)")
+    ax_e.set_xlabel("x (cm)")
+    ax_e.xaxis.set_major_formatter(mticker.ScalarFormatter(useMathText=True))
+    cur_lims = ax_e.get_ylim()
+    ax_e.set_ylim(1.e-10*cur_lims[-1], cur_lims[-1])
+
+    f.tight_layout()
+    f.savefig("profile_compare.png")
+
+
+if __name__ == "__main__":
+
+
+    data = [("res0.003km", "300 cm"),
+            ("res0.024km", "2400 cm"),
+            ("res0.192km", "0.192 km"),
+            ("res1.536km", "1.536 km"),
+            ("res12.288km", "12.288 km")]
+
+    plot_Te(data)