From 64e13044cdf3ad1b00b674f20543bce9a10e7fd0 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 29 May 2024 09:21:28 -0400 Subject: [PATCH 1/6] a setup for exploring shock burning --- Exec/science/Detonation/shock_paper/README.md | 23 ++ .../Detonation/shock_paper/det_speed_comp.py | 31 +++ .../Detonation/shock_paper/detonation.py | 105 ++++++++ .../shock_paper/inputs-det-x.subch_base | 112 ++++++++ .../Detonation/shock_paper/profile_compare.py | 60 +++++ .../Detonation/shock_paper/profiles.py | 244 ++++++++++++++++++ .../Detonation/shock_paper/show_shock_flag.py | 58 +++++ .../Detonation/shock_paper/zoom_summary.py | 74 ++++++ 8 files changed, 707 insertions(+) create mode 100644 Exec/science/Detonation/shock_paper/README.md create mode 100644 Exec/science/Detonation/shock_paper/det_speed_comp.py create mode 100644 Exec/science/Detonation/shock_paper/detonation.py create mode 100644 Exec/science/Detonation/shock_paper/inputs-det-x.subch_base create mode 100755 Exec/science/Detonation/shock_paper/profile_compare.py create mode 100755 Exec/science/Detonation/shock_paper/profiles.py create mode 100755 Exec/science/Detonation/shock_paper/show_shock_flag.py create mode 100755 Exec/science/Detonation/shock_paper/zoom_summary.py diff --git a/Exec/science/Detonation/shock_paper/README.md b/Exec/science/Detonation/shock_paper/README.md new file mode 100644 index 0000000000..650823ead6 --- /dev/null +++ b/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 | + + + + diff --git a/Exec/science/Detonation/shock_paper/det_speed_comp.py b/Exec/science/Detonation/shock_paper/det_speed_comp.py new file mode 100644 index 0000000000..d300993c9b --- /dev/null +++ b/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") diff --git a/Exec/science/Detonation/shock_paper/detonation.py b/Exec/science/Detonation/shock_paper/detonation.py new file mode 100644 index 0000000000..710371e132 --- /dev/null +++ b/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])) + diff --git a/Exec/science/Detonation/shock_paper/inputs-det-x.subch_base b/Exec/science/Detonation/shock_paper/inputs-det-x.subch_base new file mode 100644 index 0000000000..c233d2d926 --- /dev/null +++ b/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 + diff --git a/Exec/science/Detonation/shock_paper/profile_compare.py b/Exec/science/Detonation/shock_paper/profile_compare.py new file mode 100755 index 0000000000..35edb2502c --- /dev/null +++ b/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) diff --git a/Exec/science/Detonation/shock_paper/profiles.py b/Exec/science/Detonation/shock_paper/profiles.py new file mode 100755 index 0000000000..26aa70cb8f --- /dev/null +++ b/Exec/science/Detonation/shock_paper/profiles.py @@ -0,0 +1,244 @@ +#!/usr/bin/env python3 + +# Take a sequence of plotfiles and plot T and enuc vs. position + +import argparse +import re +import sys + +import matplotlib +import numpy as np +from cycler import cycler + +matplotlib.use('agg') +import math + +import matplotlib.pyplot as plt +import yt + + +## Define RGBA to HEX +def rgba_to_hex(rgba): + r = int(rgba[0]*255.0) + g = int(rgba[1]*255.0) + b = int(rgba[2]*255.0) + return f'#{r:02X}{g:02X}{b:02X}' + + +# Extract number from nuc list +def nuc_list_filter(nuc): + + match = re.search(r'\d+', nuc) + if match: + return int(match.group()) + + return 0 + + +def get_Te_profile(plotfile, plot_in_nse=False): + + ds = yt.load(plotfile, hint="castro") + + 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]) + if plot_in_nse: + in_nse = np.array(ad['in_nse'][srt]) + return time, x_coord, temp, enuc, in_nse + return time, x_coord, temp, enuc + +def get_nuc_profile(plotfile): + + ds = yt.load(plotfile, hint="castro") + + 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]) + + nuc_list = [f[1] for f in ds.field_list if f[1][0] == "X"] + nuc_list.sort(key=nuc_list_filter) + + nuc_fracs = [np.array(ad[nuc][srt]) for nuc in nuc_list] + + return time, x_coord, nuc_fracs + + +def plot_Te(prefix, nums, skip, limitlabels, xmin, xmax, plot_in_nse=False): + + f = plt.figure() + + # Get set of colors to use and apply to plot + numplots = int(len(nums) / skip) + cm = plt.get_cmap('nipy_spectral') + clist = [cm(0.95*i/numplots) for i in range(numplots + 1)] + hexclist = [rgba_to_hex(ci) for ci in clist] + + if plot_in_nse: + f.set_size_inches(7.0, 12.0) + ax_nse = f.add_subplot(311) + ax_T = f.add_subplot(312) + ax_e = f.add_subplot(313) + ax_nse.set_prop_cycle(cycler('color', hexclist)) + else: + f.set_size_inches(7.0, 9.0) + ax_T = f.add_subplot(211) + ax_e = f.add_subplot(212) + + ax_T.set_prop_cycle(cycler('color', hexclist)) + ax_e.set_prop_cycle(cycler('color', hexclist)) + + if limitlabels > 1: + skiplabels = int(numplots / limitlabels) + elif limitlabels < 0: + print("Illegal value for limitlabels: %.0f" % limitlabels) + sys.exit() + else: + skiplabels = 1 + index = 0 + + for n in range(0, len(nums), skip): + + pfile = f"{prefix}{nums[n]}" + + if plot_in_nse: + time, x, T, enuc, in_nse = get_Te_profile(pfile, plot_in_nse) + ax_nse.plot(x, in_nse) + else: + time, x, T, enuc = get_Te_profile(pfile) + + if index % skiplabels == 0: + ax_T.plot(x, T, label=f"t = {time:6.4g} s") + else: + ax_T.plot(x, T) + + ax_e.plot(x, enuc) + + index = index + 1 + + ax_T.legend(frameon=False) + ax_T.set_ylabel("T (K)") + + if xmax > 0: + ax_T.set_xlim(xmin, xmax) + ax_e.set_xlim(xmin, xmax) + if plot_in_nse: + ax_nse.set_xlim(xmin, xmax) + + ax_e.set_yscale("log") + ax_e.set_ylabel(r"$S_\mathrm{nuc}$ (erg/g/s)") + ax_e.set_xlabel("x (cm)") + + cur_lims = ax_e.get_ylim() + ax_e.set_ylim(1.e-10*cur_lims[-1], cur_lims[-1]) + + if plot_in_nse: + ax_nse.set_ylabel("IN NSE") + + f.savefig("det_Te.png") + + +def plot_nuc_frac(prefix, nums, skip, limitlabels, xmin, xmax): + + f = plt.figure() + f.set_size_inches(32.0, 20.0) + + # Get set of colors to use and apply to plot + numplots = int(len(nums) / skip) + cm = plt.get_cmap('nipy_spectral') + clist = [cm(0.95*i/numplots) for i in range(numplots + 1)] + hexclist = [rgba_to_hex(ci) for ci in clist] + + if limitlabels > 1: + skiplabels = int(numplots / limitlabels) + elif limitlabels < 0: + print("Illegal value for limitlabels: %.0f" % limitlabels) + sys.exit() + else: + skiplabels = 1 + + pfile = f"{prefix}{nums[1]}" + ds = yt.load(pfile, hint="castro") + + nuc_list = [f[1] for f in ds.field_list if f[1][0] == "X"] + nuc_list.sort(key=nuc_list_filter) + N = len(nuc_list) + + nrows = math.ceil(math.sqrt(N)) + ncols = math.ceil(math.sqrt(N)) + + for i in range(N): + ax = f.add_subplot(nrows, ncols, i+1) + ax.set_prop_cycle(cycler('color', hexclist)) + + index = 0 + for n in range(0, len(nums), skip): + + pfile = f"{prefix}{nums[n]}" + + time, x, nuc_prof = get_nuc_profile(pfile) + + if i == 0 and index % skiplabels == 0: + ax.plot(x, nuc_prof[i], label=f"t = {time:6.4g} s") + else: + ax.plot(x, nuc_prof[i]) + + index = index + 1 + + ax.legend(frameon=False) + ax.set_ylabel(nuc_list[i]) + ax.set_yscale("log") + + if xmax > 0: + ax.set_xlim(xmin, xmax) + + f.tight_layout() + f.savefig("det_nuc.png") + + +def doit(prefix, nums, skip, limitlabels, xmin, xmax, + do_nuc_fracs=False, plot_in_nse=False): + + if do_nuc_fracs: + plot_nuc_frac(prefix, nums, skip, limitlabels, xmin, xmax) + else: + plot_Te(prefix, nums, skip, limitlabels, xmin, xmax, plot_in_nse) + + +if __name__ == "__main__": + + p = argparse.ArgumentParser() + + p.add_argument("--skip", type=int, default=1, + help="interval between plotfiles") + p.add_argument("--xmin", type=float, default=0, + help="minimum x-coordinate to show") + p.add_argument("--xmax", type=float, default=-1, + help="maximum x-coordinate to show") + p.add_argument("plotfiles", type=str, nargs="+", + help="list of plotfiles to plot") + p.add_argument("--limitlabels", type=float, default=1., + help="Show all labels (default) or reduce to ~ given value") + p.add_argument("--do_nuc_fracs", dest="do_nuc_fracs", + action="store_true", + help="Plot nuc fracs, otherwise Temp and enuc plot") + p.add_argument("--plot_in_nse", dest="plot_in_nse", + action="store_true", + help="Plot in_nse quantity along with temperature and enuc") + + 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) + + doit(plot_prefix, plot_nums, args.skip, args.limitlabels, args.xmin, + args.xmax, do_nuc_fracs=args.do_nuc_fracs, plot_in_nse=args.plot_in_nse) diff --git a/Exec/science/Detonation/shock_paper/show_shock_flag.py b/Exec/science/Detonation/shock_paper/show_shock_flag.py new file mode 100755 index 0000000000..969d68d7a7 --- /dev/null +++ b/Exec/science/Detonation/shock_paper/show_shock_flag.py @@ -0,0 +1,58 @@ +#!/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(ddir): + + f = plt.figure() + + f.set_size_inches(7.5, 9.0) + + ax_T = f.add_subplot(211) + ax_e = f.add_subplot(212) + + d = detonation.Detonation(ddir) + profile = d.get_data(-1) + + ax_T.plot(profile.x, profile.T) + + ishk = profile.shock == 1.0 + ax_T.scatter(profile.x[ishk], profile.T[ishk]) + + ax_e.plot(profile.x, profile.enuc) + ax_e.scatter(profile.x[ishk], profile.enuc[ishk]) + + 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_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("shock_flag.png") + + +if __name__ == "__main__": + + #ddir = "res12.288km" + #ddir = "res1.536km" + #ddir = "res0.192km" + ddir = "res0.003km" + + plot_Te(ddir) diff --git a/Exec/science/Detonation/shock_paper/zoom_summary.py b/Exec/science/Detonation/shock_paper/zoom_summary.py new file mode 100755 index 0000000000..f0a6aa32f8 --- /dev/null +++ b/Exec/science/Detonation/shock_paper/zoom_summary.py @@ -0,0 +1,74 @@ +#!/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) + + idx = np.argmax(profile.enuc) + + xpeak = profile.x[idx] + + if n == 0: + Lx = profile.x.max() - profile.x.min() + + xmin = xpeak - 0.002 * Lx + xmax = xpeak + 0.001 * Lx + + ax_T.set_xlim(xmin, xmax) + ax_e.set_xlim(xmin, xmax) + + ax_T.scatter(profile.x, profile.T, label=label, marker="*") + ax_e.scatter(profile.x, profile.enuc, label=label, marker="*") + + 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("summary_zoom_Te.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) From 409f82e19d6d2bdccd2d4f283914cc7985b2e3e8 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 29 May 2024 09:57:46 -0400 Subject: [PATCH 2/6] add a script to setup the runs --- ....subch_base => inputs-shock-burn.template} | 9 ++-- .../Detonation/shock_paper/setup_runs.py | 48 +++++++++++++++++++ 2 files changed, 54 insertions(+), 3 deletions(-) rename Exec/science/Detonation/shock_paper/{inputs-det-x.subch_base => inputs-shock-burn.template} (91%) create mode 100755 Exec/science/Detonation/shock_paper/setup_runs.py diff --git a/Exec/science/Detonation/shock_paper/inputs-det-x.subch_base b/Exec/science/Detonation/shock_paper/inputs-shock-burn.template similarity index 91% rename from Exec/science/Detonation/shock_paper/inputs-det-x.subch_base rename to Exec/science/Detonation/shock_paper/inputs-shock-burn.template index c233d2d926..7758414c87 100644 --- a/Exec/science/Detonation/shock_paper/inputs-det-x.subch_base +++ b/Exec/science/Detonation/shock_paper/inputs-shock-burn.template @@ -7,7 +7,7 @@ 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 +amr.n_cell = @nx@ 16 16 # >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< @@ -30,6 +30,9 @@ 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 @@ -53,8 +56,8 @@ 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.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 diff --git a/Exec/science/Detonation/shock_paper/setup_runs.py b/Exec/science/Detonation/shock_paper/setup_runs.py new file mode 100755 index 0000000000..15c0b525b9 --- /dev/null +++ b/Exec/science/Detonation/shock_paper/setup_runs.py @@ -0,0 +1,48 @@ +#!/usr/bin/env python + +import os +import shutil + +# tuples of label, nx, and max_level +RES_INFO = [("12.288km", "48", "0"), + ("1.536km", "384", "0"), + ("0.192km", "3072", "0"), + ("0.024km", "6144", "1"), + ("0.003km", "12288", "2")] + +INPUTS_TEMPLATE = "inputs-shock-burn.template" + +COMMON_FILES = ["helm_table.dat", + "Castro1d.gnu.MPI.SMPLSDC.ex"] + +shock_flag = "0" +shock_thresh = "0.6666" + +def doit(): + + # read in the template + with open(INPUTS_TEMPLATE) as tf: + template = tf.readlines() + + # loop over resolutions + for label, nx, max_level in RES_INFO: + + # create output direct + odir = f"res{label}" + if shock_flag == 1: + odir += "_noshockburn" + + os.mkdir(odir) + + # copy files + for f in COMMON_FILES: + shutil.copy(f, odir) + + # modify inputs + with open(f"{odir}/inputs", "w") as fin: + for line in template: + fin.write(line.replace("@nx@", nx).replace("@nlevel@", max_level).replace("@shock_flag@", shock_flag).replace("@shock_thresh@", shock_thresh)) + + +if __name__ == "__main__": + doit() From 868b936783dee3c413f487c349c3e3e3b6473477 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 29 May 2024 10:57:04 -0400 Subject: [PATCH 3/6] more work on vis --- .../Detonation/shock_paper/det_speed_comp.py | 26 ++++++++++++++++--- .../Detonation/shock_paper/setup_runs.py | 4 +-- 2 files changed, 25 insertions(+), 5 deletions(-) mode change 100644 => 100755 Exec/science/Detonation/shock_paper/det_speed_comp.py diff --git a/Exec/science/Detonation/shock_paper/det_speed_comp.py b/Exec/science/Detonation/shock_paper/det_speed_comp.py old mode 100644 new mode 100755 index d300993c9b..fa884ef10b --- a/Exec/science/Detonation/shock_paper/det_speed_comp.py +++ b/Exec/science/Detonation/shock_paper/det_speed_comp.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python + import matplotlib.pyplot as plt import detonation @@ -5,8 +7,13 @@ runs = [("res12.288km", 12.288), ("res1.536km", 1.536), ("res0.192km", 0.192), - ("res0.024km", 0.024), - ("res0.003km", 0.003)] + ("res0.024km", 0.024)] #, +#("res0.003km", 0.003)] + +nsb_runs = [("res12.288km_noshockburn", 12.288), + ("res1.536km_noshockburn", 1.536), + ("res0.192km_noshockburn", 0.192), + ("res0.024km_noshockburn", 0.024)] #, res = [] v = [] @@ -18,14 +25,27 @@ v.append(d.v) dv.append(d.v_sigma) +nsb_res = [] +nsb_v = [] +nsb_dv = [] + +for ddir, dx in nsb_runs: + nsb_res.append(dx) + d = detonation.Detonation(ddir) + nsb_v.append(d.v) + nsb_dv.append(d.v_sigma) + fig, ax = plt.subplots() -ax.errorbar(res, v, yerr=dv, fmt="o") +ax.errorbar(res, v, yerr=dv, fmt="o", label="burning in shocks allowed") +ax.errorbar(nsb_res, nsb_v, yerr=nsb_dv, fmt="d", label="no shock burning") 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") diff --git a/Exec/science/Detonation/shock_paper/setup_runs.py b/Exec/science/Detonation/shock_paper/setup_runs.py index 15c0b525b9..42bc7ff015 100755 --- a/Exec/science/Detonation/shock_paper/setup_runs.py +++ b/Exec/science/Detonation/shock_paper/setup_runs.py @@ -15,7 +15,7 @@ COMMON_FILES = ["helm_table.dat", "Castro1d.gnu.MPI.SMPLSDC.ex"] -shock_flag = "0" +shock_flag = "1" shock_thresh = "0.6666" def doit(): @@ -29,7 +29,7 @@ def doit(): # create output direct odir = f"res{label}" - if shock_flag == 1: + if shock_flag == "1": odir += "_noshockburn" os.mkdir(odir) From 93a60429257a20a621f40176839ee3d70d8bdfdf Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 31 May 2024 15:40:19 -0400 Subject: [PATCH 4/6] a better spread --- .../Detonation/shock_paper/det_speed_comp.py | 55 +++++++++++++------ .../shock_paper/inputs-shock-burn.template | 6 +- .../Detonation/shock_paper/setup_runs.py | 16 +++--- 3 files changed, 51 insertions(+), 26 deletions(-) diff --git a/Exec/science/Detonation/shock_paper/det_speed_comp.py b/Exec/science/Detonation/shock_paper/det_speed_comp.py index fa884ef10b..7d1ae9fe9c 100755 --- a/Exec/science/Detonation/shock_paper/det_speed_comp.py +++ b/Exec/science/Detonation/shock_paper/det_speed_comp.py @@ -4,16 +4,28 @@ import detonation -runs = [("res12.288km", 12.288), - ("res1.536km", 1.536), - ("res0.192km", 0.192), - ("res0.024km", 0.024)] #, +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)] -nsb_runs = [("res12.288km_noshockburn", 12.288), - ("res1.536km_noshockburn", 1.536), - ("res0.192km_noshockburn", 0.192), - ("res0.024km_noshockburn", 0.024)] #, +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 = [] @@ -25,21 +37,32 @@ v.append(d.v) dv.append(d.v_sigma) -nsb_res = [] -nsb_v = [] -nsb_dv = [] +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 nsb_runs: - nsb_res.append(dx) +for ddir, dx in nsb1_runs: + nsb1_res.append(dx) d = detonation.Detonation(ddir) - nsb_v.append(d.v) - nsb_dv.append(d.v_sigma) + 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(nsb_res, nsb_v, yerr=nsb_dv, fmt="d", label="no shock burning") +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)") diff --git a/Exec/science/Detonation/shock_paper/inputs-shock-burn.template b/Exec/science/Detonation/shock_paper/inputs-shock-burn.template index 7758414c87..a10f50260e 100644 --- a/Exec/science/Detonation/shock_paper/inputs-shock-burn.template +++ b/Exec/science/Detonation/shock_paper/inputs-shock-burn.template @@ -1,12 +1,12 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- max_step = 200000 -stop_time = 0.02 +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 = 5.89824e7 2500 2500 +geometry.prob_hi = 1.152e8 2500 2500 amr.n_cell = @nx@ 16 16 @@ -69,7 +69,7 @@ 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.plot_per = 2.5e-3 amr.derive_plot_vars = ALL # problem initialization diff --git a/Exec/science/Detonation/shock_paper/setup_runs.py b/Exec/science/Detonation/shock_paper/setup_runs.py index 42bc7ff015..416ef4bac7 100755 --- a/Exec/science/Detonation/shock_paper/setup_runs.py +++ b/Exec/science/Detonation/shock_paper/setup_runs.py @@ -4,11 +4,13 @@ import shutil # tuples of label, nx, and max_level -RES_INFO = [("12.288km", "48", "0"), - ("1.536km", "384", "0"), - ("0.192km", "3072", "0"), - ("0.024km", "6144", "1"), - ("0.003km", "12288", "2")] +RES_INFO = [("24.0km", "48", "0"), + ("12.0km", "96", "0"), + ("6.0km", "192", "0"), + ("3.0km", "384", "0"), + ("1.5km", "768", "0"), + ("0.1875km", "6144", "0"), + ("0.0234375km", "12288", "1")] INPUTS_TEMPLATE = "inputs-shock-burn.template" @@ -16,7 +18,7 @@ "Castro1d.gnu.MPI.SMPLSDC.ex"] shock_flag = "1" -shock_thresh = "0.6666" +shock_thresh = "0.666" def doit(): @@ -30,7 +32,7 @@ def doit(): # create output direct odir = f"res{label}" if shock_flag == "1": - odir += "_noshockburn" + odir += f"_noshockburn_{shock_thresh}" os.mkdir(odir) From 0ef8d00a5c05be5d4b1d80f6e205fe398f1e6c9e Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 4 Jun 2024 14:28:45 -0400 Subject: [PATCH 5/6] more updates --- Exec/science/Detonation/shock_paper/README.md | 21 ++++++++++++------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/Exec/science/Detonation/shock_paper/README.md b/Exec/science/Detonation/shock_paper/README.md index 650823ead6..2b3d01f81b 100644 --- a/Exec/science/Detonation/shock_paper/README.md +++ b/Exec/science/Detonation/shock_paper/README.md @@ -6,18 +6,23 @@ 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: +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) | | ------------ | ----------- | ------------------- | -| 12.288 km | 48 | 1 | -| 1.536 km | 384 | 1 | -| 0.192 km | 3072 | 1 | -| 2400 cm | 6144 | 2 | -| 300 cm | 12288 | 3 | - +| 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. From 17a063380389c2756fdf4b9578235aa30da1cc1a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 4 Jun 2024 14:35:08 -0400 Subject: [PATCH 6/6] update the README --- Exec/science/Detonation/shock_paper/README.md | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/Exec/science/Detonation/shock_paper/README.md b/Exec/science/Detonation/shock_paper/README.md index 2b3d01f81b..cecf4e2236 100644 --- a/Exec/science/Detonation/shock_paper/README.md +++ b/Exec/science/Detonation/shock_paper/README.md @@ -24,5 +24,24 @@ the following resolutions into separate directories (using the 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.