Merge branch 'development' into add_subch_base_input

.clang-tidy

@@ -7,7 +7,6 @@ Checks: >
      clang-diagnostic-*,
      cppcoreguidelines-*,
        -cppcoreguidelines-avoid-c-arrays,
-       -cppcoreguidelines-avoid-goto,
        -cppcoreguidelines-avoid-magic-numbers,
        -cppcoreguidelines-avoid-non-const-global-variables,
        -cppcoreguidelines-init-variables,
@@ -17,11 +16,17 @@ Checks: >
        -cppcoreguidelines-non-private-member-variables-in-classes,
        -cppcoreguidelines-owning-memory,
        -cppcoreguidelines-pro-*,
+     misc-*,
+       -misc-const-correctness,
+       -misc-include-cleaner,
+       -misc-non-private-member-variables-in-classes,
      modernize-*,
        -modernize-avoid-c-arrays,
        -modernize-use-trailing-return-type,
        -modernize-use-using,
      performance-*,
+       -performance-avoid-endl,
+     portability-*,
      readability-*,
        -readability-avoid-const-params-in-decls,
        -readability-braces-around-statements,

.github/workflows/docs-test.yml

@@ -9,6 +9,10 @@ on:
     branches:
       - development
 
+env:
+  # enable color output from Sphinx
+  FORCE_COLOR: "1"
+
 jobs:
   docs:
     runs-on: ubuntu-latest
@@ -33,7 +37,10 @@ jobs:
       - name: Build docs
         run: |
           cd Docs/
-          make SPHINXOPTS=-v NO_DOXYGEN=TRUE html
+          # remove missing TOC entries that would otherwise be generated by
+          # doxygen to avoid warnings about missing references
+          sed -i -e 's/^   filelist$//; s/^   classlist$//' source/index.rst
+          make SPHINXOPTS='-v -W --keep-going' NO_DOXYGEN=TRUE html
 
       - name: Check links
         run: |

.github/workflows/good_defines.txt

@@ -3,6 +3,7 @@ AMREX_PARTICLES
 AMREX_SPACEDIM
 AMREX_USE_CUDA
 AMREX_USE_GPU
+AMREX_USE_HIP
 AMREX_USE_OMP
 AUX_THERMO
 BL_FORT_USE_LOWERCASE

.gitignore

@@ -72,6 +72,8 @@ Docs/source/namespacelist.rst
 Docs/source/runtime_parameters.rst
 Docs/source/*_files.rst
 Docs/source/preprocessed_files
+Docs/source/yt_example.rst
+Docs/source/yt_example_files/
 
 
 amr_diag.out

Diagnostics/DustCollapse/main.cpp

@@ -172,17 +172,17 @@ int main(int argc, char* argv[])
                 const Box& bx = mfi.tilebox();
 
 #if (AMREX_SPACEDIM == 1)
-                fdustcollapse1d(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()),
+                fdustcollapse1d(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()),
                                 BL_TO_FORTRAN_FAB(lev_data_mf[mfi]),
                                 nbins, dens.dataPtr(),
                                 imask.dataPtr(), mask_size, r1, dens_comp, &cnt);
 #elif (AMREX_SPACEDIM == 2)
 
-                fdustcollapse2d(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()),
+                fdustcollapse2d(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()),
                                 BL_TO_FORTRAN_FAB(lev_data_mf[mfi]),
                                 nbins, dens.dataPtr(), volcount.dataPtr(),
                                 imask.dataPtr(), mask_size, r1,
-                                ZFILL(level_dx), dx_fine, yctr, dens_comp);
+                                AMREX_ZFILL(level_dx), dx_fine, yctr, dens_comp);
 #else
                 fdustcollapse3d(bx.loVect(), bx.hiVect(),
                                 BL_TO_FORTRAN_FAB(lev_data_mf[mfi]),

Diagnostics/Radiation/gaussian_pulse.cpp

@@ -123,11 +123,11 @@ int main(int argc, char* argv[])
         for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) {
             const Box& bx = mfi.tilebox();
 
-            fgaussian_pulse(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()),
+            fgaussian_pulse(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()),
                             BL_TO_FORTRAN_FAB(lev_data_mf[mfi]),
                             nbins, rad_bin.dataPtr(), ncount.dataPtr(),
                             imask.dataPtr(), mask_size, r1,
-                            rad_comp, ZFILL(dx), dx_fine, xctr, yctr);
+                            rad_comp, AMREX_ZFILL(dx), dx_fine, xctr, yctr);
 
         }
 

Diagnostics/Radiation/lgt_frnt1d.cpp

@@ -114,13 +114,13 @@ int main(int argc, char* argv[])
         for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) {
             const Box& bx = mfi.tilebox();
 
-            flgt_frnt1d(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()),
+            flgt_frnt1d(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()),
                         BL_TO_FORTRAN_FAB(lev_data_mf[mfi]),
                         nbins, dens_bin.dataPtr(), vel_bin.dataPtr(),
                         pres_bin.dataPtr(), rad_bin.dataPtr(),
                         imask.dataPtr(), mask_size, r1,
                         dens_comp, xmom_comp, pres_comp, rad_comp,
-                        ZFILL(dx), dx_fine);
+                        AMREX_ZFILL(dx), dx_fine);
         }
 
         // adjust r1 for the next lowest level

Diagnostics/Radiation/rad_shock.cpp

@@ -135,7 +135,7 @@ int main(int argc, char* argv[])
         for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) {
             const Box& bx = mfi.tilebox();
 
-            fradshock(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()),
+            fradshock(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()),
                       problo.dataPtr(), probhi.dataPtr(),
                       BL_TO_FORTRAN_FAB(lev_data_mf[mfi]),
                       nbins, vars_bin.dataPtr(),

Diagnostics/Radiation/rad_source.cpp

@@ -81,7 +81,7 @@ int main(int argc, char* argv[])
 
         Real rhoe, rad;
 
-        fradsource(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()),
+        fradsource(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()),
                    BL_TO_FORTRAN_FAB(lev_data_mf[mfi]),
                    &rhoe, &rad, rhoe_comp, rad_comp);
 

Diagnostics/Radiation/rad_sphere.cpp

@@ -163,12 +163,12 @@ int main(int argc, char* argv[])
         for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) {
             const Box& bx = mfi.tilebox();
 
-            fradsphere(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()),
-                       ZFILL(problo), ZFILL(probhi),
+            fradsphere(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()),
+                       AMREX_ZFILL(problo), AMREX_ZFILL(probhi),
                        BL_TO_FORTRAN_FAB(lev_data_mf[mfi]),
                        nbins, vars_bin.dataPtr(),
                        imask.dataPtr(), mask_size, r1,
-                       ZFILL(dx), &cnt);
+                       AMREX_ZFILL(dx), &cnt);
         }
 
         // adjust r1 for the next lowest level

Diagnostics/Radiation/rhd_shocktube.cpp

@@ -150,7 +150,7 @@ int main(int argc, char* argv[])
     for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) {
         const Box& bx = mfi.tilebox();
 
-        frhdshocktube(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()),
+        frhdshocktube(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()),
                       BL_TO_FORTRAN_FAB(lev_data_mf[mfi]),
                       nbins, dens_bin.dataPtr(), vel_bin.dataPtr(),
                       pres_bin.dataPtr(), rad_bin.dataPtr(),

Docs/rp.py

@@ -7,17 +7,17 @@
 from more_itertools import unique_everseen
 
 MAIN_HEADER = """
-+--------------------------------------------+-------------------------------------------------------------+---------------+
-| parameter                                  | description                                                 | default value |
-+============================================+=============================================================+===============+
++--------------------------------------------+-------------------------------------------------------------+-----------------------------+
+| parameter                                  | description                                                 | default value               |
++============================================+=============================================================+=============================+
 """
 
 SEPARATOR = """
-+--------------------------------------------+-------------------------------------------------------------+---------------+
++--------------------------------------------+-------------------------------------------------------------+-----------------------------+
 """
 
 ENTRY = """
-| {:42} | {:59} | {:13} |
+| {:42} | {:59} | {:27} |
 """
 
 WRAP_LEN = 59

Docs/source/FlowChart.rst

@@ -309,7 +309,7 @@ In the code, the objective is to evolve the state from the old time,
    A. Create ``Sborder``, initialized from ``S_old``
 
    B. Call ``clean_state()`` to make sure the thermodynamics are in
-     sync, in particular, compute the temperature.
+      sync, in particular, compute the temperature.
 
    C. [``SHOCK_VAR``] zero out the shock flag.
 

Docs/source/Hydrodynamics.rst

@@ -342,6 +342,14 @@ accounted for in **Steps 1** and **6**. The source terms are:
    S_{{\rm ext},\rho Y_k}
    \end{array}\right)^n.
 
+.. index:: USE_SPECIES_SOURCES
+
+.. note:: To reduce memory usage, we do not include source terms for the
+   advected quantities, species, and auxiliary variables in the conserved
+   state vector by default. If your application needs external source terms for
+   these variables, set ``USE_SPECIES_SOURCES=TRUE`` when compiling so that space
+   will be allocated for them.
+
 Primitive Forms
 ===============
 
@@ -585,9 +593,6 @@ runtime parameters for hydrodynamics:
 
    See :ref:`sponge_section` for more details on the sponge.
 
--  ``castro.normalize_species``: enforce that :math:`\sum_i X_i = 1`
-   (0 or 1; default: 0)
-
 .. index:: castro.small_dens, castro.small_temp, castro.small_pres
 
 Several floors are imposed on the thermodynamic quantities to prevet unphysical

Docs/source/docutils.conf

@@ -0,0 +1,2 @@
+[parsers]
+line_length_limit = 1000000

Docs/source/faq.rst

@@ -106,10 +106,10 @@ Debugging
 
    ::
 
-           print_state(mf, IntVect(D_DECL(10, 20, 30)));
+           print_state(mf, IntVect(AMREX_D_DECL(10, 20, 30)));
 
    Here, the IntVect has the dimension that we were compiled with
-   (and this is handled through the preprocessor ``D_DECL``). In
+   (and this is handled through the preprocessor ``AMREX_D_DECL``). In
    this case, we are inspecting zone (10, 20, 30), in the global index
    space. Note that since a multifab exists only on a single level, the
    integer indices here refer to the global index space on that level.

Docs/source/index.rst

@@ -60,7 +60,7 @@ https://github.com/amrex-astro/Castro
 
    filelist
    classlist
-   .. namespacelist
+.. namespacelist
 
 .. toctree::
    :caption: References

Docs/source/io.rst

@@ -271,8 +271,6 @@ radiation quantities).
 | (where X is any of the species    | :math:`\omegadot_k = DX_k/Dt`                     |                                      |
 | defined in the network)           |                                                   |                                      |
 +-----------------------------------+---------------------------------------------------+--------------------------------------+
-| ``enuc``                          | Nuclear energy generation rate / gram             | :math:`{\rm erg~g^{-1}~s^{-1}}`      |
-+-----------------------------------+---------------------------------------------------+--------------------------------------+
 | ``rho_enuc``                      | Nuclear energy generation rate density            | :math:`{\rm erg~cm^{-3}~s^{-1}}`     |
 +-----------------------------------+---------------------------------------------------+--------------------------------------+
 | ``phiGrav``                       | Gravitational potential                           | :math:`{\rm erg~g^{-1}}`             |
@@ -393,7 +391,8 @@ Derived variables
 | ``y_velocity``,                   | :math:`\ub = (\rho \ub)/\rho`                     |                             |                                         |
 | ``z_velocity``                    |                                                   |                             |                                         |
 +-----------------------------------+---------------------------------------------------+-----------------------------+-----------------------------------------+
-
+| ``enuc``                          | Nuclear energy generation rate / gram             | ``derenuc``                 | :math:`{\rm erg~g^{-1}~s^{-1}}`         |
++-----------------------------------+---------------------------------------------------+-----------------------------+-----------------------------------------+
 
 problem-specific plotfile variables
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Exec/gravity_tests/DustCollapse/problem_bc_fill.H

@@ -30,7 +30,7 @@ void problem_bc_fill(int i, int j, int k,
     Real r = std::sqrt(x * x + y * y + z * z);
 
     // XLO
-    if (bcs(0).lo(0) == FOEXTRAP && i < domlo[0]) {
+    if (bcs(0).lo(0) == amrex::BCType::foextrap && i < domlo[0]) {
 
         int ic = domlo[0];
         Real xc = problo[0] + (static_cast<Real>(ic) + 0.5_rt) * dx[0] - problem::center[0];
@@ -50,7 +50,7 @@ void problem_bc_fill(int i, int j, int k,
     }
 
     // XHI
-    if (bcs(0).hi(0) == FOEXTRAP && i > domhi[0]) {
+    if (bcs(0).hi(0) == amrex::BCType::foextrap && i > domhi[0]) {
 
         int ic = domhi[0];
         Real xc = problo[0] + (static_cast<Real>(ic) + 0.5_rt) * dx[0] - problem::center[0];
@@ -71,7 +71,7 @@ void problem_bc_fill(int i, int j, int k,
 
 #if AMREX_SPACEDIM >= 2
     // YLO
-    if (bcs(0).lo(1) == FOEXTRAP && j < domlo[1]) {
+    if (bcs(0).lo(1) == amrex::BCType::foextrap && j < domlo[1]) {
 
         int jc = domlo[1];
         Real yc = problo[1] + (static_cast<Real>(jc) + 0.5_rt) * dx[1] - problem::center[1];
@@ -91,7 +91,7 @@ void problem_bc_fill(int i, int j, int k,
     }
 
     // YHI
-    if (bcs(0).hi(1) == FOEXTRAP && j > domhi[1]) {
+    if (bcs(0).hi(1) == amrex::BCType::foextrap && j > domhi[1]) {
 
         int jc = domhi[1];
         Real yc = problo[1] + (static_cast<Real>(jc) + 0.5_rt) * dx[1] - problem::center[1];
@@ -112,7 +112,7 @@ void problem_bc_fill(int i, int j, int k,
 
 #if AMREX_SPACEDIM == 3
     // ZLO
-    if (bcs(0).lo(2) == FOEXTRAP && k < domlo[2]) {
+    if (bcs(0).lo(2) == amrex::BCType::foextrap && k < domlo[2]) {
 
         int kc = domlo[2];
         Real zc = problo[2] + (static_cast<Real>(kc) + 0.5_rt) * dx[2] - problem::center[2];
@@ -132,7 +132,7 @@ void problem_bc_fill(int i, int j, int k,
     }
 
     // ZHI
-    if (bcs(0).hi(2) == FOEXTRAP && k > domhi[2]) {
+    if (bcs(0).hi(2) == amrex::BCType::foextrap && k > domhi[2]) {
 
         int kc = domhi[2];
         Real zc = problo[2] + (static_cast<Real>(kc) + 0.5_rt) * dx[2] - problem::center[2];

Exec/gravity_tests/hydrostatic_adjust/problem_bc_fill.H

@@ -14,7 +14,7 @@ void problem_bc_fill(int i, int j, int k,
 
     if (AMREX_SPACEDIM == 1) {
 
-        if (bcs(UMX).hi(0) == FOEXTRAP && i > domhi[0]) {
+        if (bcs(UMX).hi(0) == amrex::BCType::foextrap && i > domhi[0]) {
 
             Real vel = amrex::max(state(i,j,k,UMX) / state(i,j,k,URHO), 0.0);
             state(i,j,k,URHO)  = problem::hse_rho_top;
@@ -36,7 +36,7 @@ void problem_bc_fill(int i, int j, int k,
     }
     else if (AMREX_SPACEDIM == 2) {
 
-        if (bcs(UMY).hi(1) == FOEXTRAP && j > domhi[1]) {
+        if (bcs(UMY).hi(1) == amrex::BCType::foextrap && j > domhi[1]) {
 
             Real vel = amrex::max(state(i,j,k,UMY) / state(i,j,k,URHO), 0.0);
             state(i,j,k,URHO)  = problem::hse_rho_top;
@@ -58,7 +58,7 @@ void problem_bc_fill(int i, int j, int k,
     }
     else {
 
-        if (bcs(UMZ).hi(2) == FOEXTRAP && k > domhi[2]) {
+        if (bcs(UMZ).hi(2) == amrex::BCType::foextrap && k > domhi[2]) {
 
             Real vel = amrex::max(state(i,j,k,UMZ) / state(i,j,k,URHO), 0.0);
             state(i,j,k,URHO)  = problem::hse_rho_top;

Exec/hydro_tests/double_mach_reflection/problem_bc_fill.H

@@ -18,7 +18,7 @@ void problem_bc_fill(int i, int j, int k,
     // x boundaries
     //-------------------------------------------------------------------------
 
-    if ((bcs(URHO).lo(0) == EXT_DIR || bcs(URHO).lo(0) == FOEXTRAP) && i < domlo[0]) {
+    if ((bcs(URHO).lo(0) == amrex::BCType::ext_dir || bcs(URHO).lo(0) == amrex::BCType::foextrap) && i < domlo[0]) {
         state(i,j,k,URHO) = problem::rho_l;
         state(i,j,k,UMX) = problem::rho_l * problem::u_l;
         state(i,j,k,UMY) = problem::rho_l * problem::v_l;
@@ -43,7 +43,7 @@ void problem_bc_fill(int i, int j, int k,
     Real y = problo[1] + dx[1] * (static_cast<Real>(j) + 0.5e0_rt);
 
     // YLO
-    if ((bcs(URHO).lo(1) == EXT_DIR || bcs(URHO).lo(1) == FOEXTRAP) && j < domlo[1]) {
+    if ((bcs(URHO).lo(1) == amrex::BCType::ext_dir || bcs(URHO).lo(1) == amrex::BCType::foextrap) && j < domlo[1]) {
         if (x < 1.e0_rt / 6.e0_rt) {
             // ICs
             state(i,j,k,URHO) = problem::rho_l;
@@ -71,7 +71,7 @@ void problem_bc_fill(int i, int j, int k,
     }
 
     // YHI
-    if ((bcs(URHO).hi(2) == EXT_DIR || bcs(URHO).hi(2) == FOEXTRAP) && j > domhi[1]) {
+    if ((bcs(URHO).hi(2) == amrex::BCType::ext_dir || bcs(URHO).hi(2) == amrex::BCType::foextrap) && j > domhi[1]) {
         state(i,j,k,URHO ) = 0.e0_rt;
         state(i,j,k,UMX  ) = 0.e0_rt;
         state(i,j,k,UMY  ) = 0.e0_rt;

Exec/hydro_tests/oddeven/inputs.2d

@@ -22,6 +22,7 @@ castro.do_hydro = 1
 castro.do_react = 0
 castro.ppm_type = 0
 castro.hybrid_riemann = 1
+castro.shock_detection_threshold = 0.5
 
 # TIME STEP CONTROL