Merge branch 'development' into wdmerger_hybrid

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
+