a setup for exploring shock burning (#2857)

Exec/science/Detonation/shock_paper/README.md

@@ -0,0 +1,47 @@
+# 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 the script `setup_runs.py` will setup a suite of simulations with
+the following resolutions into separate directories (using the
+`inputs-shock-burn.template`):
+
+
+| resolution   |  base grid  |  levels (4x jumps)  |
+| ------------ | ----------- | ------------------- |
+|       24 km  |       48    |         1           |
+|       12 km  |       96    |         1           |
+|        6 km  |      192    |         1           |
+|        3 km  |      384    |         1           |
+|      1.5 km  |      768    |         1           |
+|   0.1875 km  |     6144    |         1           |
+|  2343.74 cm  |    12288    |         2           |
+
+you can set the value of the shock detection threshold there
+and the directory names will reflect that setting.
+
+## plotting
+
+The following scripts can make useful plots (some use the
+`detonation.py` module as support):
+
+* `det_speed_comp.py` : plot detonation speed vs. resolution, using
+  simple differencing to estimate the detonation speed from the last 3
+  plotfiles.
+
+* `profile_compare.py` : given a list of pruns (from different
+  resolutions), make a plot showing all of their profiles together.
+
+* `profiles.py` : for a single run, plot profiles of T and enuc for
+  several times.
+
+* `show_shock_flag.py` : simply plot the shock variable on top of T
+  and enuc profiles.
+
+* `zoom_summary.py` : given a list of runs (from different
+  resolutions), plot the last plotfile from each run, zoomed in on
+  where the peak energy generation is.

Exec/science/Detonation/shock_paper/det_speed_comp.py

@@ -0,0 +1,74 @@
+#!/usr/bin/env python
+
+import matplotlib.pyplot as plt
+
+import detonation
+
+runs = [("res24.0km", 24),
+        ("res12.0km", 12),
+        ("res6.0km", 6),
+        ("res3.0km", 3),
+        ("res1.5km", 1.5),
+        ("res0.1875km", 0.1875)] #,
+       #("res0.024km", 0.024)] #,
+#("res0.003km", 0.003)]
+
+nsb1_runs = [("res24.0km_noshockburn_1", 24),
+             ("res12.0km_noshockburn_1", 12),
+             ("res6.0km_noshockburn_1", 6),
+             ("res3.0km_noshockburn_1", 3),
+             ("res1.5km_noshockburn_1", 1.5),
+             ("res0.1875km_noshockburn_1", 0.1875)] #,
+
+nsb23_runs = [("res24.0km_noshockburn_0.666", 24),
+              ("res12.0km_noshockburn_0.666", 12),
+              ("res6.0km_noshockburn_0.666", 6),
+              ("res3.0km_noshockburn_0.666", 3),
+              ("res1.5km_noshockburn_0.666", 1.5),
+              ("res0.1875km_noshockburn_0.666", 0.1875)] #,
+
+res = []
+v = []
+dv = []
+
+for ddir, dx in runs:
+    res.append(dx)
+    d = detonation.Detonation(ddir)
+    v.append(d.v)
+    dv.append(d.v_sigma)
+
+nsb23_res = []
+nsb23_v = []
+nsb23_dv = []
+
+for ddir, dx in nsb23_runs:
+    nsb23_res.append(dx)
+    d = detonation.Detonation(ddir)
+    nsb23_v.append(d.v)
+    nsb23_dv.append(d.v_sigma)
+
+nsb1_res = []
+nsb1_v = []
+nsb1_dv = []
+
+for ddir, dx in nsb1_runs:
+    nsb1_res.append(dx)
+    d = detonation.Detonation(ddir)
+    nsb1_v.append(d.v)
+    nsb1_dv.append(d.v_sigma)
+
+
+fig, ax = plt.subplots()
+
+ax.errorbar(res, v, yerr=dv, fmt="o", label="burning in shocks allowed")
+ax.errorbar(nsb23_res, nsb23_v, yerr=nsb23_dv, fmt="d", label="no shock burning (f=2/3)")
+ax.errorbar(nsb1_res, nsb1_v, yerr=nsb1_dv, fmt="d", label="no shock burning (f=1)")
+
+ax.set_xlabel(r"$\Delta x$ (km)")
+ax.set_ylabel(r"$v_\mathrm{det}$ (cm / s)")
+
+ax.legend()
+
+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-shock-burn.template

@@ -0,0 +1,115 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 200000
+stop_time =  0.03
+
+# 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     = 1.152e8  2500  2500
+amr.n_cell           = @nx@    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.disable_shock_burning = @shock_flag@
+castro.shock_detection_threshold = @shock_thresh@
+
+# 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       = @nlevel@       # maximum level number allowed
+amr.ref_ratio       = 4 4 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.5e-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
+