Merge branch 'development' into massive_star_screening

.github/workflows/dependencies_hip.sh

@@ -1,17 +1,26 @@
 #!/usr/bin/env bash
 #
-# Copyright 2020-2022 The AMReX Community
+# Copyright 2020 The AMReX Community
 #
 # License: BSD-3-Clause-LBNL
 # Authors: Axel Huebl
 
+# search recursive inside a folder if a file contains tabs
+#
+# @result 0 if no files are found, else 1
+#
+
 set -eu -o pipefail
 
+# `man apt.conf`:
+#   Number of retries to perform. If this is non-zero APT will retry
+#   failed files the given number of times.
+echo 'Acquire::Retries "3";' | sudo tee /etc/apt/apt.conf.d/80-retries
 
 # Ref.: https://rocmdocs.amd.com/en/latest/Installation_Guide/Installation-Guide.html#ubuntu
 curl -O https://repo.radeon.com/rocm/rocm.gpg.key
 sudo apt-key add rocm.gpg.key
-echo 'deb [arch=amd64] https://repo.radeon.com/rocm/apt/debian/ ubuntu main' \
+echo "deb [arch=amd64] https://repo.radeon.com/rocm/apt/${1-debian}/ ubuntu main" \
   | sudo tee /etc/apt/sources.list.d/rocm.list
 echo 'export PATH=/opt/rocm/llvm/bin:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin:$PATH' \
   | sudo tee -a /etc/profile.d/rocm.sh
@@ -34,12 +43,14 @@ sudo apt-get install -y --no-install-recommends \
     roctracer-dev   \
     rocprofiler-dev \
     rocrand-dev     \
-    rocprim-dev
+    rocprim-dev     \
+    hiprand-dev
 
 # activate
 #
 source /etc/profile.d/rocm.sh
 hipcc --version
+hipconfig --full
 which clang
 which clang++
 which flang

CHANGES.md

@@ -1,3 +1,7 @@
+# 23.12
+
+  * The radiation solver port to C++ has been completed (#2638, #2648)
+
 # 23.11
 
   * Problem GNUmakefiles have been standardized and now allow for the

Docs/source/io.rst

@@ -562,7 +562,15 @@ By default, 4 output files are created:
 
        The species masses are given in units of solar masses.
 
-Some problems have custom versions of the diagnostics with additional information.
+``Castro/Util/scripts/diag_parser.py`` contains Python code for parsing
+these output files into Numpy arrays.  Usage instructions are included
+in the file, along with an example script at
+``Castro/Util/scripts/plot_species.py``.  This reads a
+``species_diag.out`` file provided on the command line and makes a plot
+of the total mass fractions over time.
+
+Some problems have custom versions of the diagnostics with additional
+information.  These are not currently supported by the Python parser.
 
 
 .. _sec:parallel_io:

Exec/Make.Castro

@@ -270,6 +270,10 @@ ifeq ($(USE_ROTATION), TRUE)
   DEFINES += -DROTATION
 endif
 
+ifeq ($(USE_SPECIES_SOURCES), TRUE)
+  DEFINES += -DCONS_SPECIES_HAVE_SOURCES
+endif
+
 ifeq ($(USE_PARTICLES), TRUE)
   Bdirs += Source/particles
 endif

Exec/reacting_tests/nse_test/problem_initialize_state_data.H

@@ -5,6 +5,7 @@
 #include <eos.H>
 #include <network.H>
 #if defined(NSE_TABLE)
+#include <nse_table_type.H>
 #include <nse_table.H>
 #elif defined(NSE_NET)
 #include <nse_solver.H>
@@ -64,14 +65,15 @@ void problem_initialize_state_data (int i, int j, int k,
     burn_state.y_e = ye;
 
 #if defined(NSE_TABLE)
-    Real abar;
-    Real dq;
-    Real dyedt;
-    Real dabardt;
-    Real dbeadt;
-    Real e_nu;
-
-    nse_interp(T, problem::rho0, ye, abar, dq, dyedt, dabardt, dbeadt, e_nu, burn_state.xn);
+    nse_table_t nse_state;
+    nse_state.T = T;
+    nse_state.rho = problem::rho0;
+    nse_state.Ye = ye;
+    nse_interp(nse_state);
+
+    for (int n = 0; n < NumSpec; ++n) {
+        burn_state.xn[n] = nse_state.X[n];
+    }
 #elif defined(NSE_NET)
     Real eps = 1.e-10_rt;
     bool input_ye_is_valid = true;
@@ -105,16 +107,16 @@ void problem_initialize_state_data (int i, int j, int k,
     for (int n = 0; n < NumSpec; n++) {
         burn_state.xn[n] /= sumX;
     }
-    
+
 #ifdef NSE_NET
     state(i,j,k,UMUP) = burn_state.mu_p;
     state(i,j,k,UMUN) = burn_state.mu_n;
 #endif
-    
+
 #ifdef AUX_THERMO
     burn_state.aux[AuxZero::iye] = ye;
-    burn_state.aux[AuxZero::iabar] = abar;
-    burn_state.aux[AuxZero::ibea] = dq;
+    burn_state.aux[AuxZero::iabar] = nse_state.abar;
+    burn_state.aux[AuxZero::ibea] = nse_state.bea;
 #endif
 
     eos(eos_input_rt, burn_state);

Exec/science/massive_star/_prob_params

@@ -1,7 +1,10 @@
 model_name      character       ""        y
 
+interpolate_pres   int          0         y
+
 perturb_model   int             0         y
 velpert_amplitude  real         1.e5      y
 velpert_scale   real            1.e7      y
 
 isi28           int            -1
+

Exec/science/massive_star/analysis/massive_star_multi.py

@@ -1,22 +1,24 @@
 #!/usr/bin/env python3
 
+import argparse
+import os
+
 import matplotlib
-matplotlib.use('agg')
+import matplotlib.pyplot as plt
 
-import os
-import sys
 import yt
-import matplotlib.pyplot as plt
-import numpy as np
-from functools import reduce
+from yt.frontends.boxlib.api import CastroDataset
+# assume that our data is in CGS
+from yt.units import cm
+
+matplotlib.use('agg')
 
-from mpl_toolkits.axes_grid1 import ImageGrid
 
-# assume that our data is in CGS
-from yt.units import cm, amu
-from yt.frontends.boxlib.api import CastroDataset
 
-def make_plot(plotfile, prefix="plot", size=(19.2, 10.8), cbar_location="right"):
+
+def make_plot(plotfile, prefix="plot", width_frac=0.1,
+              layout=(1, 4),
+              size=(19.2, 10.8), cbar_location="right"):
 
     ds = CastroDataset(plotfile)
 
@@ -36,8 +38,6 @@ def make_plot(plotfile, prefix="plot", size=(19.2, 10.8), cbar_location="right")
     fig = plt.figure()
     fig.set_size_inches(size)
 
-    width_frac = 0.1
-
     fields = ["MachNumber", "magvort", "abar", "enuc"]
 
     sp = yt.SlicePlot(ds, "theta", fields,
@@ -75,17 +75,39 @@ def make_plot(plotfile, prefix="plot", size=(19.2, 10.8), cbar_location="right")
 
     sp.set_axes_unit("cm")
 
-    fig = sp.export_to_mpl_figure((1, len(fields)), cbar_location=cbar_location, cbar_pad="5%")
+    fig = sp.export_to_mpl_figure((layout[0], layout[1]), axes_pad=(1.0, 0.4),
+                                  cbar_location=cbar_location, cbar_pad="2%")
 
     fig.subplots_adjust(left=0.05, right=0.95, bottom=0.025, top=0.975)
-    fig.text(0.02, 0.02, "time = {:8.5f} s".format(float(ds.current_time)), transform=fig.transFigure)
+    fig.text(0.02, 0.02, f"time = {float(ds.current_time):8.2f} s",
+             transform=fig.transFigure)
     fig.set_size_inches(size)
+    fig.tight_layout()
+    extra = ""
+    if layout[0] >= layout[1]:
+        extra = "_vertical"
 
-    fig.savefig(f"{prefix}_{os.path.basename(plotfile)}_slice.png", pad_inches=0)
+    fig.savefig(f"{prefix}_{os.path.basename(plotfile)}_w{width_frac:04.2f}{extra}.pdf", pad_inches=0)
 
 
 if __name__ == "__main__":
 
-    plotfile = sys.argv[1]
-
-    make_plot(plotfile, "all")
+    p = argparse.ArgumentParser()
+    p.add_argument("--vertical", action="store_true",
+                   help="plot 2x2 or 1x4")
+    p.add_argument("--width_fraction", type=float, default=0.1,
+                   help="fraction of domain to show")
+    p.add_argument("plotfile", type=str, nargs=1,
+                   help="plotfile to plot")
+
+    args = p.parse_args()
+    plotfile = args.plotfile[0]
+
+    if args.vertical:
+        size = (7.5, 11.0)
+        layout = (2, 2)
+    else:
+        size = (19.2, 8.5)
+        layout = (1, 4)
+
+    make_plot(plotfile, "all", layout=layout, width_frac=args.width_fraction, size=size)

Exec/science/massive_star/problem_initialize_state_data.H

@@ -5,6 +5,7 @@
 #include <eos.H>
 #include <model_parser.H>
 #ifdef NSE_TABLE
+#include <nse_table_type.H>
 #include <nse_table_check.H>
 #endif
 
@@ -33,6 +34,7 @@ void problem_initialize_state_data (int i, int j, int k,
 
     state(i,j,k,URHO) = interpolate(dist, model::idens);
     state(i,j,k,UTEMP) = interpolate(dist, model::itemp);
+    Real pres = interpolate(dist, model::ipres);
     for (int n = 0; n < NumSpec; n++) {
         state(i,j,k,UFS+n) = amrex::max(interpolate(dist, model::ispec+n), small_x);
     }
@@ -79,21 +81,16 @@ void problem_initialize_state_data (int i, int j, int k,
 
     if (nse_check) {
 
-        Real abar;
-        Real dq;
-        Real dyedt;
-        Real dabardt;
-        Real dbeadt;
-        Real e_nu;
-        Real xn[NumSpec];
+        nse_table_t nse_state;
+        nse_state.T = state(i,j,k,UTEMP);
+        nse_state.rho = state(i,j,k,URHO);
+        nse_state.Ye = state(i,j,k,UFX+AuxZero::iye);
+        nse_interp(nse_state);
 
-        nse_interp(state(i,j,k,UTEMP), state(i,j,k,URHO),
-                   state(i,j,k,UFX+AuxZero::iye), abar, dq, dyedt, dabardt, dbeadt, e_nu, xn);
-
-        state(i,j,k,UFX+AuxZero::iabar) = abar;
-        state(i,j,k,UFX+AuxZero::ibea) = dq;
+        state(i,j,k,UFX+AuxZero::iabar) = nse_state.abar;
+        state(i,j,k,UFX+AuxZero::ibea) = nse_state.bea;
         for (int n = 0; n < NumSpec; n++) {
-            state(i,j,k,UFS+n) = amrex::max(xn[n], small_x);
+            state(i,j,k,UFS+n) = amrex::max(nse_state.X[n], small_x);
         }
 
         // we just got new X's, so we need to re-renormalize them
@@ -116,22 +113,29 @@ void problem_initialize_state_data (int i, int j, int k,
         }
     }
 
-    // now call the EOS -- reload the eos_state since composition may
-    // have changed
+    // now call the EOS -- we need to use an eos_t here since we are
+    // going to take pressure as an input
 
-    burn_state.rho = state(i,j,k,URHO);
-    burn_state.T = state(i,j,k,UTEMP);
+    eos_t eos_state;
+    eos_state.rho = state(i,j,k,URHO);
+    eos_state.T = state(i,j,k,UTEMP);
+    eos_state.p = pres;
     for (int n = 0; n < NumSpec; n++) {
-        burn_state.xn[n] = state(i,j,k,UFS+n);
+        eos_state.xn[n] = state(i,j,k,UFS+n);
     }
     for (int n = 0; n < NumAux; n++) {
-        burn_state.aux[n] = state(i,j,k,UFX+n);
+        eos_state.aux[n] = state(i,j,k,UFX+n);
     }
 
-    eos(eos_input_rt, burn_state);
+    if (problem::interpolate_pres == 1) {
+        eos(eos_input_rp, eos_state);
+        state(i,j,k,UTEMP) = eos_state.T;
+    } else {
+        eos(eos_input_rt, eos_state);
+    }
 
-    state(i,j,k,UEINT) = state(i,j,k,URHO) * burn_state.e;
-    state(i,j,k,UEDEN) = state(i,j,k,URHO) * burn_state.e;
+    state(i,j,k,UEINT) = state(i,j,k,URHO) * eos_state.e;
+    state(i,j,k,UEDEN) = state(i,j,k,URHO) * eos_state.e;
 
     // add density scaling to composition
 

Source/driver/Castro.cpp

@@ -1260,8 +1260,7 @@ Castro::initData ()
              {
                const Box& box = mfi.validbox();
 
-               tmp.resize(box, 1);
-               Elixir elix_tmp = tmp.elixir();
+               tmp.resize(box, 1, The_Async_Arena());
                auto tmp_arr = tmp.array();
 
                make_fourth_in_place(box, Sborder.array(mfi), tmp_arr, domain_lo, domain_hi);
@@ -1287,8 +1286,7 @@ Castro::initData ()
            {
              const Box& box = mfi.growntilebox(2);
 
-             tmp.resize(box, 1);
-             Elixir elix_tmp = tmp.elixir();
+             tmp.resize(box, 1, The_Async_Arena());
              auto tmp_arr = tmp.array();
 
              make_cell_center_in_place(box, Sborder.array(mfi), tmp_arr, domain_lo, domain_hi);
@@ -1336,8 +1334,7 @@ Castro::initData ()
            {
              const Box& box = mfi.validbox();
 
-             tmp.resize(box, 1);
-             Elixir elix_tmp = tmp.elixir();
+             tmp.resize(box, 1, The_Async_Arena());
              auto tmp_arr = tmp.array();
 
              make_fourth_in_place(box, Sborder.array(mfi), tmp_arr, domain_lo, domain_hi);
@@ -3950,8 +3947,7 @@ Castro::computeTemp(
       compute_lap_term(bx0, Stemp.array(mfi), Eint_lap.array(mfi), UEINT,
                        domain_lo, domain_hi);
 
-      tmp.resize(bx, 1);
-      Elixir elix_tmp = tmp.elixir();
+      tmp.resize(bx, 1, The_Async_Arena());
       auto tmp_arr = tmp.array();
 
       make_cell_center_in_place(bx, Stemp.array(mfi), tmp_arr, domain_lo, domain_hi);
@@ -4067,8 +4063,7 @@ Castro::computeTemp(
 
       const Box& bx = mfi.tilebox();
 
-      tmp.resize(bx, 1);
-      Elixir elix_tmp = tmp.elixir();
+      tmp.resize(bx, 1, The_Async_Arena());
       auto tmp_arr = tmp.array();
 
       // only temperature