Merge branch 'development' into add_cubic_interp

.clang-tidy

@@ -13,7 +13,6 @@ Checks: >
        -cppcoreguidelines-init-variables,
        -cppcoreguidelines-interfaces-global-init,
        -cppcoreguidelines-macro-usage,
-       -cppcoreguidelines-no-malloc,
        -cppcoreguidelines-non-private-member-variables-in-classes,
        -cppcoreguidelines-owning-memory,
        -cppcoreguidelines-pro-*,
@@ -25,7 +24,6 @@ Checks: >
      modernize-*,
        -modernize-avoid-c-arrays,
        -modernize-use-trailing-return-type,
-       -modernize-use-using,
      performance-*,
        -performance-avoid-endl,
      portability-*,

.github/workflows/check_pr_branch.yml

@@ -0,0 +1,20 @@
+name: check PR branch
+
+on:
+  pull_request:
+    types:
+      - opened
+      - synchronize
+      - reopened
+      - edited
+
+jobs:
+  check-PR-branch:
+    runs-on: ubuntu-latest
+    steps:
+      - name: PRs should not target main
+        run: |
+          if [[ "${{ github.base_ref }}" == "main" ]]; then
+            echo 'Pull requests must not be made against main. Please target development instead.'
+            exit 1
+          fi

Exec/Make.Castro

@@ -357,7 +357,7 @@ include $(AMREX_HOME)/Tools/GNUMake/Make.rules
 clean::
 	$(SILENT) $(RM) extern.F90 prob_params_auto.F90 extern_parameters.H extern_parameters_F.H extern_parameters.cpp
 	$(SILENT) $(RM) AMReX_buildInfo.cpp
-	$(SILENT) $(RM) $(CASTRO_AUTO_SOURCE_DIR)/*.H $(CASTRO_AUTO_SOURCE_DIR)/*.[fF]90
+	$(SILENT) $(RM) $(CASTRO_AUTO_SOURCE_DIR)/*.H $(CASTRO_AUTO_SOURCE_DIR)/*.[fF]90 $(CASTRO_AUTO_SOURCE_DIR)/*.cpp
 
 # these files are now created directly in the CASTRO_AUTO_SOURCE_DIR so eventually we
 # can get rid of explicitly removing them (this is for backwards compatibility)

Exec/science/flame_wave/inputs_3d.testsuite.gpu

@@ -100,7 +100,7 @@ problem.T_hi = 2.e8
 problem.T_lo   = 8.e6
 
 problem.H_star = 2000.e0
-problem.atm_delta  = 50.0
+problem.atm_delta  = 100.0
 
 problem.fuel1_name = "helium-4"
 problem.fuel1_frac = 1.0e0

Exec/science/massive_star/analysis/massive_star_slice_3d.py

@@ -0,0 +1,107 @@
+#!/usr/bin/env python3
+
+import os
+import sys
+
+import matplotlib
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import ImageGrid
+
+import yt
+from yt.frontends.boxlib.api import CastroDataset
+
+matplotlib.use('agg')
+
+
+plotfile = sys.argv[1]
+ds = CastroDataset(plotfile)
+
+t_drive = 0.0
+if "[*] castro.drive_initial_convection_tmax" in ds.parameters:
+    t_drive = ds.parameters["[*] castro.drive_initial_convection_tmax"]
+elif "castro.drive_initial_convection_tmax" in ds.parameters:
+    t_drive = ds.parameters["castro.drive_initial_convection_tmax"]
+print(t_drive)
+
+domain_frac = 0.2
+
+xmin = ds.domain_left_edge[0]
+xmax = ds.domain_right_edge[0]
+xctr = 0.5 * (xmin + xmax)
+L_x = xmax - xmin
+xmin = xctr - 0.5 * domain_frac * L_x
+xmax = xctr + 0.5 * domain_frac * L_x
+L_x = xmax - xmin
+
+ymin = ds.domain_left_edge[1]
+ymax = ds.domain_right_edge[1]
+yctr = 0.5 * (ymin + ymax)
+L_y = ymax - ymin
+ymin = yctr - 0.5 * domain_frac * L_y
+ymax = yctr + 0.5 * domain_frac * L_y
+L_y = ymax - ymin
+
+zmin = ds.domain_left_edge[2]
+zmax = ds.domain_right_edge[2]
+zctr = 0.5 * (zmin + zmax)
+L_z = zmax - zmin
+zmin = zctr - 0.5 * domain_frac * L_z
+zmax = zctr + 0.5 * domain_frac * L_z
+L_z = zmax - zmin
+
+fig = plt.figure()
+
+
+fields = ["MachNumber", "magvort", "abar", "enuc"]
+
+grid = ImageGrid(fig, 111, nrows_ncols=(2, (len(fields)+1)//2),
+                 axes_pad=0.75, cbar_pad=0.05, label_mode="L", cbar_mode="each")
+
+
+for i, f in enumerate(fields):
+
+    sp = yt.SlicePlot(ds, "y", f,
+                      center=[xctr, yctr, zctr], width=[L_x, L_z],
+                      fontsize="12")
+
+    sp.set_buff_size((2400,2400))
+    sp.swap_axes()
+
+    if f == "Ye":
+        sp.set_zlim(f, 0.46, 0.5)
+        sp.set_log(f, False)
+        sp.set_cmap(f, "magma_r")
+    elif f == "abar":
+        sp.set_log(f, False)
+        sp.set_cmap(f, "viridis")
+    elif f == "enuc":
+        sp.set_log(f, True, linthresh=1.e12)
+        sp.set_zlim(f, -1.e20, 1.e20)
+        sp.set_cmap(f, "bwr")
+    elif f == "MachNumber":
+        sp.set_zlim(f, 1.e-4, 0.3)
+        sp.set_cmap(f, "plasma")
+    elif f == "magvel":
+        sp.set_zlim(f, 100.0, 2.e7)
+        sp.set_cmap(f, "viridis")
+    elif f == "magvort":
+        sp.set_cmap(f, "magma")
+        sp.set_zlim(f, 1.e-2, 5)
+
+    sp.set_axes_unit("km")
+
+    plot = sp.plots[f]
+    plot.figure = fig
+    plot.axes = grid[i].axes
+    plot.cax = grid.cbar_axes[i]
+    if i < len(fields)-1:
+        grid[i].axes.xaxis.offsetText.set_visible(False)
+
+    sp._setup_plots()
+
+fig.text(0.02, 0.02, f"$t - \\tau_\\mathrm{{drive}}$ = {float(ds.current_time) - t_drive:6.1f} s",
+         transform=fig.transFigure)
+
+fig.set_size_inches(11, 10)
+plt.tight_layout()
+plt.savefig(f"{os.path.basename(plotfile)}_slice.png")

Exec/science/massive_star/inputs_3d.nse

@@ -3,18 +3,18 @@
 amr.plot_files_output = 1
 amr.checkpoint_files_output = 1
 
-max_step = 50
+max_step = 5000000
 stop_time = 360000
 
-geometry.is_periodic = 0 0
+geometry.is_periodic = 0 0 0
 geometry.coord_sys = 0         # r-z coordinates
 
 geometry.prob_lo   =  0.    0.    0.
 geometry.prob_hi   =  1.6384e10 1.6384e10 1.6384e10
 
 amr.n_cell         = 512 512 512
 
-amr.max_level      = 3      # maximum level number allowed
+amr.max_level      = 2      # maximum level number allowed
 
 castro.lo_bc       =  2 2 2
 castro.hi_bc       =  2 2 2
@@ -50,21 +50,21 @@ castro.sponge_upper_density = 1.e3
 castro.sponge_lower_density = 1.e2
 castro.sponge_timescale     = 1.e-3
 
-castro.cfl            = 0.4     # cfl number for hyperbolic system
-castro.init_shrink    = 0.01     # scale back initial timestep by this factor
+castro.cfl            = 0.5     # cfl number for hyperbolic system
+castro.init_shrink    = 0.1     # scale back initial timestep by this factor
 castro.change_max     = 1.2    # factor by which dt is allowed to change each timestep
 castro.sum_interval   = 10       # timesteps between computing and printing volume averages
 
 #castro.dtnuc_e = 0.25
 #castro.dtnuc_X = 0.25
 
-amr.ref_ratio       = 4 2 2 2 # refinement ratio
+amr.ref_ratio       = 4 4 2 2 2 # refinement ratio
 amr.regrid_int      = 10000   # how often to regrid
-amr.n_error_buf     = 4 2 2 2 # number of buffer cells in error est
+amr.n_error_buf     = 4 4 2 2 2 # number of buffer cells in error est
 amr.grid_eff        = 0.7     # what constitutes an efficient grid
 
 amr.check_file      = massive_star_chk     # root name of checkpoint file
-amr.check_int       = 200     # number of timesteps between checkpoints
+amr.check_int       = 10     # number of timesteps between checkpoints
 
 amr.plot_file       = massive_star_plt     # root name of plot file
 amr.plot_per = 5.0
@@ -76,10 +76,6 @@ amr.small_plot_per = 0.5
 amr.small_plot_vars = density Temp in_nse
 amr.derive_small_plot_vars = abar Ye enuc MachNumber magvel magvort
 
-#amr.checkpoint_files_output = 0
-#amr.plot_files_output = 0
-#amr.smallplot_files_output = 0
-
 fab.format = NATIVE_32
 
 castro.plot_per_is_exact = 0
@@ -88,6 +84,8 @@ castro.plot_per_is_exact = 0
 amr.max_grid_size   = 64       # maximum grid size allowed -- used to control parallelism
 amr.blocking_factor = 32       # block factor in grid generation
 
+amr.subcycling_mode = None
+
 amr.v               = 1       # control verbosity in Amr.cpp
 castro.v            = 1       # control verbosity in Castro.cpp
 
@@ -116,17 +114,13 @@ castro.drive_initial_convection_tmax = 50
 
 # refinement
 
-amr.refinement_indicators = denerr denerr2 denerr3
+amr.refinement_indicators = denerr denerr3
 
 amr.refine.denerr.max_level = 1
 amr.refine.denerr.value_greater = 2.e3
 amr.refine.denerr.field_name = density
 
-amr.refine.denerr2.max_level = 2
-amr.refine.denerr2.value_greater = 3.e4
-amr.refine.denerr2.field_name = density
-
-amr.refine.denerr3.max_level = 3
+amr.refine.denerr3.max_level = 2
 amr.refine.denerr3.value_greater = 3.e5
 amr.refine.denerr3.field_name = density
 
@@ -146,7 +140,15 @@ network.Si_nse = 0.02
 network.C_nse = 1.0
 network.O_nse = 1.0
 
-integrator.ode_max_steps = 500000
+integrator.ode_max_steps = 5000
+
+integrator.use_burn_retry = 1
+integrator.retry_swap_jacobian = 1
+
+integrator.retry_rtol_spec = 1.e-5
+integrator.retry_atol_spec = 1.e-5
+integrator.retry_rtol_enuc = 1.e-5
+integrator.retry_atol_enuc = 1.e-5
 
 network.small_x = 1.e-10
 

Source/driver/Castro_io.cpp

@@ -7,7 +7,6 @@
 #include <iostream>
 #include <string>
 #include <ctime>
-#include <filesystem>
 
 #include <AMReX_Utility.H>
 #include <Castro.H>
@@ -579,7 +578,7 @@ Castro::writeJobInfo (const std::string& dir, const Real io_time)
   jobInfoFile << "output date / time: "
               << std::put_time(std::localtime(&now), "%c\n") << "\n";
 
-  jobInfoFile << "output dir:         " << std::filesystem::current_path() << "\n";
+  jobInfoFile << "output dir:         " << amrex::FileSystem::CurrentPath() << "\n";
 
   jobInfoFile << "I/O time (s):       " << io_time << "\n";
 

Source/gravity/Gravity.cpp

@@ -2966,17 +2966,21 @@ Gravity::set_mass_offset (Real time, bool multi_level) const
         {
             for (int lev = 0; lev <= parent->finestLevel(); lev++) {
                 auto* cs = dynamic_cast<Castro*>(&parent->getLevel(lev));
-        if (cs != nullptr) {
-            mass_offset += cs->volWgtSum("density", time);
-        } else {
+                if (cs != nullptr) {
+                    mass_offset += cs->volWgtSum("density", time);
+                } else {
                     amrex::Abort("unable to access volWgtSum");
                 }
             }
         }
         else
         {
             auto* cs = dynamic_cast<Castro*>(&parent->getLevel(0));
-            mass_offset = cs->volWgtSum("density", time, false, false);  // do not mask off fine grids
+            if (cs != nullptr) {
+                mass_offset = cs->volWgtSum("density", time, false, false);  // do not mask off fine grids
+            } else {
+                amrex::Abort("unable to access volWgtSum");
+            }
         }
 
         mass_offset = mass_offset / geom.ProbSize();

Source/radiation/Radiation.cpp

@@ -5,8 +5,6 @@
 #include <rad_util.H>
 #include <filt_prim.H>
 
-#include <AMReX_PROB_AMR_F.H>
-
 #include <opacity.H>
 
 #include <iostream>

Source/reactions/Castro_react.cpp

@@ -533,7 +533,9 @@ Castro::react_state(Real time, Real dt)
 #endif
     int num_failed = 0;
 
-    // why no omp here?
+#ifdef _OPENMP
+#pragma omp parallel reduction(+:num_failed)
+#endif
     for (MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi)
     {
         const Box& bx = mfi.growntilebox(ng);

Source/sdc/Castro_sdc.cpp

@@ -1,5 +1,4 @@
 #include <Castro.H>
-// #include <Castro_hydro_F.H>
 #include <Castro_sdc_util.H>
 
 using namespace amrex;

Source/sources/Castro_sources.cpp

@@ -529,6 +529,15 @@ Castro::pre_advance_operators (Real time, Real dt)  // NOLINT(readability-conver
 
     advance_status status {};
 
+    // If we are using gravity, solve for the potential and
+    // gravitational field.  note: since reactions don't change
+    // density, we can do this before or after the burn.
+
+#ifdef GRAVITY
+    construct_old_gravity(time);
+#endif
+
+
     // If we are Strang splitting the reactions, do the old-time contribution now.
 
 #ifndef TRUE_SDC
@@ -541,11 +550,6 @@ Castro::pre_advance_operators (Real time, Real dt)  // NOLINT(readability-conver
 #endif
 #endif
 
-    // If we are using gravity, solve for the potential and gravitational field.
-
-#ifdef GRAVITY
-    construct_old_gravity(time);
-#endif
 
     // Initialize the new-time data. This copy needs to come after all Strang-split operators.