Updating nova setup (#3031)

This updates the nova problem to match the A4 model in
the nova paper that is about to be submitted.

Exec/science/nova/GNUmakefile

@@ -1,7 +1,7 @@
 PRECISION        = DOUBLE
 PROFILE          = FALSE
 DEBUG            = FALSE
-DIM              = 2
+DIM              = 3
 
 COMP	         = gnu
 
@@ -10,7 +10,7 @@ USE_OMP          = FALSE
 USE_CUDA         = FALSE
 
 USE_REACT        = TRUE
-USE_DIFFUSION    = TRUE
+USE_DIFFUSION    = FALSE
 USE_GRAV         = TRUE
 
 CASTRO_HOME := ../../..
@@ -22,7 +22,7 @@ MAX_NPTS_MODEL = 20000
 EOS_DIR     := helmholtz
 
 # This sets the network directory in $(MICROPHYSICS_HOME)/networks
-NETWORK_DIR := nova2
+NETWORK_DIR := nova
 
 # Thi sets the conductivity EOS directory in $(MICROPHYSICS_HOME)/conductivity
 CONDUCTIVITY_DIR := stellar

Exec/science/nova/_prob_params

@@ -9,6 +9,10 @@ velpert_amplitude     real            1.e2_rt       y
 
 num_vortices          integer         1             y
 
+num_vortices_x        integer         1             y
+
+num_vortices_y        integer         1             y
+
 max_num_vortices      integer         10            n
 
 xloc_vortices         real            0.0_rt        n     10
@@ -19,7 +23,9 @@ y_h1                 real             0.0_rt        n
 
 width                real            0.0_rt         y
 
-L                    real            0.0_rt         n
+L_x                  real            0.0_rt         n
+
+L_y                  real            0.0_rt         n
 
 amplitude            real            0.0_rt         y
 

Exec/science/nova/analysis/cei_energy.py

@@ -0,0 +1,172 @@
+import yt
+import os
+import re
+import numpy as np
+from yt.utilities.parallel_tools.parallel_analysis_interface import communication_system
+import matplotlib.pyplot as plt
+
+yt.enable_parallelism()
+
+plt.rcParams['font.size'] = 21
+plt.rcParams['font.family'] = 'serif'
+plt.rcParams["axes.labelsize"] = 21
+plt.rcParams['mathtext.fontset'] = 'stix'
+plt.rcParams['savefig.bbox']='tight'
+
+curr_dir = os.getcwd()
+plt_dir_a04 = [f.name for f in os.scandir(curr_dir) if (f.is_dir() and re.search("^plt[0-9]{5,8}$",f.name))]
+plt_dir_a04.sort(key=lambda x: int(x[3:]))
+
+dir_paths_a04 = []
+for dir in plt_dir_a04:
+    dir_paths_a04.append(os.path.join(curr_dir, dir))
+
+index_a04 = dict(zip(plt_dir_a04, [i for i in range(len(plt_dir_a04))]))
+
+curr_dir = "../nova_t7_0.8km_ext_pert"
+plt_dir_a08 = [f.name for f in os.scandir(curr_dir) if (f.is_dir() and re.search("^plt[0-9]{5,8}$",f.name))]
+plt_dir_a08.sort(key=lambda x: int(x[3:]))
+
+dir_paths_a08 = []
+for dir in plt_dir_a08:
+    dir_paths_a08.append(os.path.join(curr_dir, dir))
+
+index_a08 = dict(zip(plt_dir_a08, [i for i in range(len(plt_dir_a08))]))
+
+curr_dir = "../nova_t7_0.2km"
+plt_dir_b02 = [f.name for f in os.scandir(curr_dir) if (f.is_dir() and re.search("^plt[0-9]{5,8}$",f.name))]
+plt_dir_b02.sort(key=lambda x: int(x[3:]))
+
+dir_paths_b02 = []
+for dir in plt_dir_b02:
+    dir_paths_b02.append(os.path.join(curr_dir, dir))
+
+index_b02 = dict(zip(plt_dir_b02, [i for i in range(len(plt_dir_b02))]))
+
+curr_dir = "../nova_t7_0.4km"
+plt_dir_b04 = [f.name for f in os.scandir(curr_dir) if (f.is_dir() and re.search("^plt[0-9]{5,8}$",f.name))]
+plt_dir_b04.sort(key=lambda x: int(x[3:]))
+
+dir_paths_b04 = []
+for dir in plt_dir_b04:
+    dir_paths_b04.append(os.path.join(curr_dir, dir))
+
+index_b04 = dict(zip(plt_dir_b04, [i for i in range(len(plt_dir_b04))]))
+
+ts_a08 = yt.DatasetSeries(dir_paths_a08)
+ts_a04 = yt.DatasetSeries(dir_paths_a04)
+
+ts_b04 = yt.DatasetSeries(dir_paths_b04[0:138])
+ts_b02 = yt.DatasetSeries(dir_paths_b02[0:267])
+
+N_a08 = len(plt_dir_a08)
+time_a08 = np.zeros(N_a08)
+qvalue_a08 = np.zeros(N_a08)
+
+N_a04 = len(plt_dir_a04)
+time_a04 = np.zeros(N_a04)
+qvalue_a04 = np.zeros(N_a04)
+
+N_b04 = len(plt_dir_b04[0:138])
+time_b04 = np.zeros(N_b04)
+qvalue_b04 = np.zeros(N_b04)
+
+N_b02 = len(plt_dir_b02[0:267])
+time_b02 = np.zeros(N_b02)
+qvalue_b02 = np.zeros(N_b02)
+
+comm = communication_system.communicators[-1]
+
+for ds in ts_a04.piter():
+
+    id_i = index_a04[ds.basename]
+
+    sp = ds.all_data()
+
+    max_enuc = sp.max("enuc")
+    time_a04[id_i] = ds.current_time.value
+    qvalue_a04[id_i] = max_enuc
+
+    ds.index.clear_all_data()
+
+comm.barrier()
+
+for ds in ts_a08.piter():
+
+    id_i = index_a08[ds.basename]
+
+    sp = ds.all_data()
+
+    max_enuc = sp.max("enuc")
+    time_a08[id_i] = ds.current_time.value
+    qvalue_a08[id_i] = max_enuc
+
+    ds.index.clear_all_data()
+
+for ds in ts_b02.piter():
+
+    id_i = index_b02[ds.basename]
+
+    sp = ds.all_data()
+
+    max_enuc = sp.max("enuc")
+    time_b02[id_i] = ds.current_time.value
+    qvalue_b02[id_i] = max_enuc
+
+    ds.index.clear_all_data()
+
+for ds in ts_b04.piter():
+
+    id_i = index_b04[ds.basename]
+
+    sp = ds.all_data()
+
+    max_enuc = sp.max("enuc")
+    time_b04[id_i] = ds.current_time.value
+    qvalue_b04[id_i] = max_enuc
+
+    ds.index.clear_all_data()
+
+time_a04[:] = comm.mpi_allreduce(time_a04[:], op="sum")
+qvalue_a04[:] = comm.mpi_allreduce(qvalue_a04[:], op="sum")
+
+time_a08[:] = comm.mpi_allreduce(time_a08[:], op="sum")
+qvalue_a08[:] = comm.mpi_allreduce(qvalue_a08[:], op="sum")
+
+time_b02[:] = comm.mpi_allreduce(time_b02[:], op="sum")
+qvalue_b02[:] = comm.mpi_allreduce(qvalue_b02[:], op="sum")
+
+time_b04[:] = comm.mpi_allreduce(time_b04[:], op="sum")
+qvalue_b04[:] = comm.mpi_allreduce(qvalue_b04[:], op="sum")
+
+if yt.is_root():
+
+    per = 10
+    t_a04 = time_a04[::per]
+    t_a08 = time_a08[::per]
+    t_b02 = time_b02[::per]
+    t_b04 = time_b04[::per]
+
+    q_a04 = qvalue_a04[::per]
+    q_a08 = qvalue_a08[::per]
+    q_b02 = qvalue_b02[::per]
+    q_b04 = qvalue_b04[::per]
+
+    fig = plt.figure()
+    ax = fig.add_subplot(111)
+    ax.plot(t_a04, q_a04, label="Model A4", linestyle='--')
+    ax.plot(t_a08, q_a08, label="Model A8", linestyle='--')
+    ax.plot(t_b02, q_b02, label="Model B2", linestyle='-', linewidth=3.0)
+    ax.plot(t_b04, q_b04, label="Model B4", linestyle='-', linewidth=3.0)
+
+    ax.set_xlabel(r'$t\,[s]$')
+    ax.set_ylabel(r'$\dot{e}_{\mathrm{nuc},\mathrm{max}}\,[\mathrm{erg}\cdot \mathrm{g}^{-1}\,\mathrm{s}^{-1}]$')
+    ax.set_yscale('log')
+    ax.set_ylim(1.0e13, 1.0e19)
+    ax.legend(loc='upper left')
+
+    fig.set_size_inches(8, 6)
+    ax.tick_params(labelsize=21)
+    plt.tight_layout()
+
+    fig.savefig('q_series_04_ext.pdf', bbox_inches="tight")

Exec/science/nova/analysis/diag_parser.py

@@ -0,0 +1,99 @@
+"""Helper functions for working with Castro diagnostic files (*_diag.out)
+
+To use these in a standalone script, you can do one of the following:
+
+* append $CASTRO_HOME/Util/scripts to sys.path at the top of your script:
+    sys.path.append("<path to Castro>/Util/scripts")
+
+* add a symlink to this file in the same directory as your script:
+    $ ln -s "$CASTRO_HOME/Util/scripts/diag_parser.py" .
+
+* copy this file into the same directory as your script
+
+Then you can do `from diag_parser import deduplicate, read_diag_file`.
+"""
+
+from pathlib import Path
+
+import numpy as np
+
+""" Format notes
+files are opened in Castro.cpp, data is written in sum_integrated_quantities.cpp
+
+data_logs[0]: grid_diag.out
+intwidth, fixwidth, datwidth*
+
+data_logs[1]: gravity_diag.out
+- this was previously missing the last column number (8), which we handle for
+  backwards compatibility
+intwidth, fixwidth, datwidth, datwidth, datwidth, datwidth, datwidth, datwidth
+
+data_logs[2]: species_diag.out
+intwidth, fixwidth, datwidth*
+
+data_logs[3]: amr_diag.out
+- if compiled with GPU support, this will have two additional integer fields at
+  the end with size `datwidth` for the GPU memory usage
+- column 5 (max number of subcycles) is an integer
+intwidth, fixwidth, fixwidth, intwidth, fixwidth, datwidth
+"""
+
+datwidth = 25  # Floating point data in scientific notation
+fixwidth = 25  # Floating point data not in scientific notation
+intwidth = 12  # Integer data
+
+# Any additional columns after these are assumed to be floating point values in
+# scientific notation (amr_diag.out gets special handling)
+FIELD_WIDTHS = {
+    "grid_diag.out": [intwidth, fixwidth],
+    "gravity_diag.out": [intwidth, fixwidth] + [datwidth] * 6,
+    "species_diag.out": [intwidth, fixwidth],
+    "amr_diag.out": [intwidth, fixwidth, fixwidth, intwidth, fixwidth, datwidth],
+}
+
+
+def read_diag_file(file_path):
+    """Reads a Castro diagnostic file into a numpy structured array.
+
+    Currently only supports the default files that Castro generates.
+    """
+    if not isinstance(file_path, Path):
+        file_path = Path(file_path)
+    filetype = file_path.name
+    if filetype not in FIELD_WIDTHS:
+        raise ValueError("Unsupported file name")
+    widths = FIELD_WIDTHS[filetype]
+    with open(file_path, "r") as f:
+        # try getting the number of columns from the first line
+        first_line = f.readline().rstrip("\n")
+        if filetype == "gravity_diag.out":
+            # gravity_diag.out is missing the last column number, but it
+            # fortunately has a fixed number of columns
+            num_columns = 8
+        else:
+            num_columns = int(first_line.split()[-1])
+        # pad out the widths list on the right if necessary
+        widths.extend([datwidth] * (num_columns - len(widths)))
+        # infer datatypes from the widths
+        dtypes = [int if w == intwidth else float for w in widths]
+        # amr_diag.out has several integer columns with long names
+        if filetype == "amr_diag.out":
+            dtypes[4] = int  # max number of subcycles
+            if num_columns >= 8:
+                dtypes[6] = int  # maximum gpu memory used
+                dtypes[7] = int  # minimum gpu memory free
+        # already read the first header line, so we don't need to skip any rows
+        data = np.genfromtxt(
+            f, delimiter=widths, comments="#", dtype=dtypes, names=True
+        )
+    return data
+
+
+def deduplicate(data):
+    """Deduplicate based on the timestep, keeping the only last occurrence."""
+    # get the unique indices into the reversed timestep array, so we find the
+    # final occurrence of each timestep
+    _, rev_indices = np.unique(data["TIMESTEP"][::-1], return_index=True)
+    # np.unique() sorts by value, so we don't need to un-reverse rev_indices
+    unique_indices = data.shape[0] - rev_indices - 1
+    return data[unique_indices]

Exec/science/nova/analysis/lateral_dens_04_ext.py

@@ -0,0 +1,60 @@
+import yt
+import os
+import re
+import numpy as np
+from yt.utilities.parallel_tools.parallel_analysis_interface import communication_system
+import matplotlib.pyplot as plt
+from yt.frontends.boxlib.api import CastroDataset
+
+plt.rcParams['font.size'] = 12
+plt.rcParams['font.family'] = 'serif'
+plt.rcParams['mathtext.fontset'] = 'stix'
+plt.rcParams['savefig.bbox']='tight'
+
+curr_dir = os.getcwd()
+
+plt_dir = [f.name for f in os.scandir(curr_dir) if (f.is_dir() and re.search("^plt+[0-9]{5,8}$",f.name))]
+plt_dir.sort(key=lambda x: int(x[3:]))
+
+dir_paths = []
+for dir in plt_dir:
+    dir_paths.append(os.path.join(curr_dir, dir))
+
+fig = plt.figure()
+ax = fig.add_subplot(111)
+ax.set_xlabel(r'$y\,[\mathrm{km}]$')
+ax.set_ylabel(r'$\rho\,[\mathrm{g/cm^3}]$')
+
+times = [300, 931, 1524, 1532]
+
+for time in times:
+
+    if time == 1532:
+        ds =  ds = CastroDataset(dir_paths[-2])
+    else:
+        id = int(time/0.5)
+        ds = CastroDataset(dir_paths[id])
+
+    sp = ds.all_data()
+    res = 960
+    profile = yt.create_profile(
+        sp,
+        [("boxlib", "y")], [("boxlib", "density")],
+        weight_field=('boxlib', 'cell_volume'),
+        accumulation=False,
+        n_bins=res,
+        )
+
+    y = profile.x.to("km")
+    val = profile[("boxlib", "density")]
+    ax.plot(y, val, label="${:.2f}\,[s]$".format(ds.current_time.value))
+
+ax.legend()
+
+fig.set_size_inches(5, 7)
+ax.tick_params(labelsize=14)
+ax.set_ylim(1.0e-1, 5.0e5)
+ax.set_yscale('log')
+plt.tight_layout()
+
+fig.savefig('lateral_dens_04_ext.pdf', bbox_inches="tight")