From 25cc6d8f8594a613c88682da052e1147c45bd64b Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Fri, 8 Mar 2024 10:26:48 -0500 Subject: [PATCH] update detonation problem for nse_net (#2765) --- Exec/science/Detonation/GNUmakefile | 2 +- Exec/science/Detonation/analysis/profiles.py | 48 +++++++++++++------ Exec/science/Detonation/inputs-det-x.nse_net | 8 ++-- .../problem_initialize_state_data.H | 10 ++-- 4 files changed, 44 insertions(+), 24 deletions(-) diff --git a/Exec/science/Detonation/GNUmakefile b/Exec/science/Detonation/GNUmakefile index 256b2c173b..2783b93351 100644 --- a/Exec/science/Detonation/GNUmakefile +++ b/Exec/science/Detonation/GNUmakefile @@ -18,7 +18,7 @@ EOS_DIR := helmholtz # This sets the EOS directory in $(MICROPHYSICS_HOME)/networks ifeq ($(USE_NSE_NET), TRUE) - NETWORK_DIR := ase + NETWORK_DIR := subch_base SCREEN_METHOD := chabrier1998 else NETWORK_DIR := aprox19 diff --git a/Exec/science/Detonation/analysis/profiles.py b/Exec/science/Detonation/analysis/profiles.py index 2fd0f6ae3a..26aa70cb8f 100755 --- a/Exec/science/Detonation/analysis/profiles.py +++ b/Exec/science/Detonation/analysis/profiles.py @@ -35,7 +35,7 @@ def nuc_list_filter(nuc): return 0 -def get_Te_profile(plotfile): +def get_Te_profile(plotfile, plot_in_nse=False): ds = yt.load(plotfile, hint="castro") @@ -48,10 +48,11 @@ def get_Te_profile(plotfile): 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") @@ -72,19 +73,27 @@ def get_nuc_profile(plotfile): return time, x_coord, nuc_fracs -def plot_Te(prefix, nums, skip, limitlabels, xmin, xmax): +def plot_Te(prefix, nums, skip, limitlabels, xmin, xmax, plot_in_nse=False): f = plt.figure() - f.set_size_inches(7.0, 9.0) - - ax_T = f.add_subplot(211) - ax_e = f.add_subplot(212) # 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)) @@ -101,7 +110,11 @@ def plot_Te(prefix, nums, skip, limitlabels, xmin, xmax): pfile = f"{prefix}{nums[n]}" - time, x, T, enuc = get_Te_profile(pfile) + 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") @@ -117,6 +130,9 @@ def plot_Te(prefix, nums, skip, limitlabels, xmin, xmax): 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)") @@ -125,8 +141,8 @@ def plot_Te(prefix, nums, skip, limitlabels, xmin, xmax): cur_lims = ax_e.get_ylim() ax_e.set_ylim(1.e-10*cur_lims[-1], cur_lims[-1]) - if xmax > 0: - ax_e.set_xlim(xmin, xmax) + if plot_in_nse: + ax_nse.set_ylabel("IN NSE") f.savefig("det_Te.png") @@ -189,12 +205,13 @@ def plot_nuc_frac(prefix, nums, skip, limitlabels, xmin, xmax): f.savefig("det_nuc.png") -def doit(prefix, nums, skip, limitlabels, xmin, xmax, do_nuc_fracs=False): +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_Te(prefix, nums, skip, limitlabels, xmin, xmax, plot_in_nse) if __name__ == "__main__": @@ -214,6 +231,9 @@ def doit(prefix, nums, skip, limitlabels, xmin, xmax, do_nuc_fracs=False): 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() @@ -221,4 +241,4 @@ def doit(prefix, nums, skip, limitlabels, xmin, xmax, do_nuc_fracs=False): 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) + args.xmax, do_nuc_fracs=args.do_nuc_fracs, plot_in_nse=args.plot_in_nse) diff --git a/Exec/science/Detonation/inputs-det-x.nse_net b/Exec/science/Detonation/inputs-det-x.nse_net index cab9513b75..d67a09d289 100644 --- a/Exec/science/Detonation/inputs-det-x.nse_net +++ b/Exec/science/Detonation/inputs-det-x.nse_net @@ -32,7 +32,7 @@ castro.small_temp = 1.e7 castro.riemann_solver = 0 -castro.time_integration_method = 0 +castro.time_integration_method = 3 # TIME STEP CONTROL castro.cfl = 0.5 # cfl number for hyperbolic system @@ -104,10 +104,10 @@ integrator.rtol_spec = 1.e-6 integrator.atol_spec = 1.e-6 integrator.rtol_enuc = 1.e-6 -integrator.jacobian = 1 +integrator.jacobian = 2 integrator.ode_max_steps = 1500000 -nse.nse_molar_independent = 1 +nse.nse_molar_independent = 0 nse.nse_dx_independent = 1 nse.ase_tol = 0.1 -nse.nse_skip_molar = 1 \ No newline at end of file +nse.nse_skip_molar = 0 \ No newline at end of file diff --git a/Exec/science/Detonation/problem_initialize_state_data.H b/Exec/science/Detonation/problem_initialize_state_data.H index 3b08ce86a9..f594e8f15e 100644 --- a/Exec/science/Detonation/problem_initialize_state_data.H +++ b/Exec/science/Detonation/problem_initialize_state_data.H @@ -62,12 +62,12 @@ void problem_initialize_state_data (int i, int j, int k, burn_state.mu_p = problem::mu_p; burn_state.mu_n = problem::mu_n; - if (burn_state.T > 2.0e9) { - auto nse_state = get_actual_nse_state(burn_state); + if (burn_state.T > 2.5e9) { + auto nse_state = get_actual_nse_state(burn_state, 1.0e-10_rt, true); - for (int n = 0; n < NumSpec; ++n) { - state(i,j,k,UFS+n) = state(i,j,k,URHO) * nse_state.xn[n]; - } + // for (int n = 0; n < NumSpec; ++n) { + // state(i,j,k,UFS+n) = state(i,j,k,URHO) * nse_state.xn[n]; + // } state(i,j,k,UMUP) = nse_state.mu_p; state(i,j,k,UMUN) = nse_state.mu_n;