From 25afc9bfe07cfce29e383fcc6a16730e1e8460f3 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 26 Jul 2024 11:38:54 -0400 Subject: [PATCH 1/3] a setup for exploring shock burning (#2857) --- Exec/science/Detonation/shock_paper/README.md | 47 ++++ .../Detonation/shock_paper/det_speed_comp.py | 74 ++++++ .../Detonation/shock_paper/detonation.py | 105 ++++++++ .../shock_paper/inputs-shock-burn.template | 115 +++++++++ .../Detonation/shock_paper/profile_compare.py | 60 +++++ .../Detonation/shock_paper/profiles.py | 244 ++++++++++++++++++ .../Detonation/shock_paper/setup_runs.py | 50 ++++ .../Detonation/shock_paper/show_shock_flag.py | 58 +++++ .../Detonation/shock_paper/zoom_summary.py | 74 ++++++ 9 files changed, 827 insertions(+) create mode 100644 Exec/science/Detonation/shock_paper/README.md create mode 100755 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-shock-burn.template 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/setup_runs.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..cecf4e2236 --- /dev/null +++ b/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. 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 100755 index 0000000000..7d1ae9fe9c --- /dev/null +++ b/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") 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-shock-burn.template b/Exec/science/Detonation/shock_paper/inputs-shock-burn.template new file mode 100644 index 0000000000..a10f50260e --- /dev/null +++ b/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 + 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/setup_runs.py b/Exec/science/Detonation/shock_paper/setup_runs.py new file mode 100755 index 0000000000..416ef4bac7 --- /dev/null +++ b/Exec/science/Detonation/shock_paper/setup_runs.py @@ -0,0 +1,50 @@ +#!/usr/bin/env python + +import os +import shutil + +# tuples of label, nx, and max_level +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" + +COMMON_FILES = ["helm_table.dat", + "Castro1d.gnu.MPI.SMPLSDC.ex"] + +shock_flag = "1" +shock_thresh = "0.666" + +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 += f"_noshockburn_{shock_thresh}" + + 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() 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 0f6eead3a57588546f61e68d64c5e1e1fbe52729 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 26 Jul 2024 12:47:10 -0400 Subject: [PATCH 2/3] fix the subch_planar diag for cylindrical coords (#2885) the case of projecting against div{U} was not including area factors --- Exec/science/subch_planar/Problem_Derive.cpp | 21 ++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/Exec/science/subch_planar/Problem_Derive.cpp b/Exec/science/subch_planar/Problem_Derive.cpp index 2e97187e44..3275166238 100644 --- a/Exec/science/subch_planar/Problem_Derive.cpp +++ b/Exec/science/subch_planar/Problem_Derive.cpp @@ -273,12 +273,15 @@ void ca_dergradpoverp1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nc Real vm = dat(i,j-1,k,UMY) / dat(i,j-1,k,URHO); Real v0 = dat(i,j,k,UMY) / dat(i,j,k,URHO); + Real du_x{}; + Real dv_y{}; + // construct div{U} if (coord_type == 0) { // Cartesian - div_u += 0.5_rt * (up - um) * dxinv; - div_u += 0.5_rt * (vp - vm) * dyinv; + du_x = 0.5_rt * (up - um) * dxinv; + dv_y = 0.5_rt * (vp - vm) * dyinv; } else if (coord_type == 1) { @@ -287,8 +290,8 @@ void ca_dergradpoverp1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nc Real rm = (i - 1 + 0.5_rt) * dx[0]; Real rp = (i + 1 + 0.5_rt) * dx[0]; - div_u += 0.5_rt * (rp * up - rm * um) / (rc * dx[0]) + - 0.5_rt * (vp - vm) * dyinv; + du_x = 0.5_rt * (rp * up - rm * um) / (rc * dx[0]); + dv_y = 0.5_rt * (vp - vm) * dyinv; #ifndef AMREX_USE_GPU } else { @@ -296,6 +299,8 @@ void ca_dergradpoverp1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nc #endif } + div_u = du_x + dv_y; + // we need to compute p in the full stencil Real p_ip1{}; @@ -399,12 +404,12 @@ void ca_dergradpoverp1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nc Real dP_y = 0.5_rt * (p_jp1 - p_jm1); //Real gradPdx_over_P = std::sqrt(dP_x * dP_x + dP_y * dP_y + dP_z * dP_z) / dat(i,j,k,QPRES); - Real du_x = std::min(up - um, 0.0); - Real dv_y = std::min(vp - vm, 0.0); + Real cdu_x = std::min(du_x, 0.0); + Real cdv_y = std::min(dv_y, 0.0); - Real divu_mag = std::sqrt(du_x * du_x + dv_y * dv_y + 1.e-30); + Real divu_mag = std::sqrt(cdu_x * cdu_x + cdv_y * cdv_y + 1.e-30); - Real gradPdx_over_P = std::abs(dP_x * du_x + dP_y * dv_y) / divu_mag; + Real gradPdx_over_P = std::abs(dP_x * cdu_x + dP_y * cdv_y) / divu_mag; gradPdx_over_P /= p_zone; der(i,j,k,0) = gradPdx_over_P; From c8ef18138b904b9bbb8eb0246463edc121f3088e Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 26 Jul 2024 19:13:05 -0400 Subject: [PATCH 3/3] fix true-SDC compilation due to namespace changes in Microphysics (#2932) --- Source/sdc/sdc_react_util.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/sdc/sdc_react_util.H b/Source/sdc/sdc_react_util.H index dffd1dd1f0..b174dd8d03 100644 --- a/Source/sdc/sdc_react_util.H +++ b/Source/sdc/sdc_react_util.H @@ -35,7 +35,7 @@ single_zone_react_source(burn_t& burn_state, eos(eos_input_re, burn_state); // eos_get_small_temp(&small_temp); - burn_state.T = amrex::min(MAX_TEMP, amrex::max(burn_state.T, small_temp)); + burn_state.T = std::clamp(burn_state.T, small_temp, integrator_rp::MAX_TEMP); Array1D ydot;