AMReX-Astro · BenWibking · Oct 19, 2025 · Oct 19, 2025 · Oct 19, 2025 · Oct 19, 2025
diff --git a/.github/workflows/good_defines.txt b/.github/workflows/good_defines.txt
@@ -24,6 +24,7 @@ SCREENING
 SCREEN_METHOD
 SDC
 SIMPLIFIED_SDC
+SINGLE_PRECISION_JACOBIAN
 STRANG
 TRUE_SDC
 _OPENMP

diff --git a/AGENTS.md b/AGENTS.md
@@ -0,0 +1,19 @@
+# Repository Guidelines
+
+## Project Structure & Module Organization
+The repository is organized around modular physics kernels in `EOS/`, `networks/`, `integration/`, and companion directories such as `conductivity/`, `neutrinos/`, `opacity/`, and `rates/`. Shared math utilities and build scripts remain in `util/`, while reusable constants sit in `constants/`. Interfaces consumed by AMReX-based application codes are in `interfaces/`. Unit regression assets reside under `unit_test/`, and Sphinx documentation sources are in `Docs/` with release history tracked in `CHANGES.md`.
+
+## Build, Test, and Development Commands
+Export `AMREX_HOME` and `MICROPHYSICS_HOME` before building so the GNU make stubs locate AMReX sources. Build a standalone unit test in place, e.g. `cd unit_test/burn_cell && make -j4`, which yields `main3d.gnu.ex`. Select physics at compile time with flags like `make NETWORK_DIR=aprox19 EOS_DIR=helmholtz`, and reset a directory via `make clean`.
+
+## Coding Style & Naming Conventions
+C++ sources follow `.editorconfig`: four-space indentation, LF line endings, and tabs only inside Makefiles. Headers keep the uppercase `.H` suffix and implementations use `.cpp`. Favor the existing snake_case routine pattern (`nse_check`, `burn_cell`) and guard autogenerated directories such as `util/autodiff`. Run `.clang-tidy` with `make USE_CLANG_TIDY=TRUE` on significant changes.
+
+## Testing Guidelines
+Each directory beneath `unit_test/` ships a `GNUmakefile`, `inputs_*` controls, and README notes; run the built binary with the matching inputs file (for example `./main3d.gnu.ex inputs_aprox13`). For new coverage, clone an existing `test_*` layout, describe the scenario, and rerun the most relevant case before posting. Capture scratch artifacts via `.gitignore` rather than committing them.
+
+## Commit & Pull Request Guidelines
+Develop on topic branches and open pull requests against `development`; merges to `main` occur monthly. Keep commit subjects concise and imperative, mirroring history such as `add NSE network compatibility docs (#1852)`. Update `CHANGES.md` for user-facing fixes, link issues, and note required AMReX revisions. PRs should summarize physics choices, list tests run, and attach plots or logs when behavior shifts.
+
+## Documentation & Automation
+Sphinx sources in `Docs/` back the published guide; update them with code changes and ensure `make -C Docs html` finishes cleanly. GitHub Actions publishes the result and exercises unit tests, so keep workflows green by limiting optional dependencies and recording new Python requirements in `requirements.txt`.
diff --git a/Docs/source/ode_integrators.rst b/Docs/source/ode_integrators.rst
@@ -113,6 +113,23 @@ For the other networks (usually pynucastro networks), the implementation is
 provided in ``Microphysics/util/linpack.H`` and is templated on the number
 of equations.  Pivoting can be disabled by setting ``integrator.linalg_do_pivoting=0``.
 
+.. index:: USE_SINGLE_PRECISION_JACOBIAN
+
+.. tip::
+
+   The storage for the Jacobian can take up the most memory when
+   integrating the reaction system.  It is possible to store the
+   Jacobian as single-precision, by building with:
+
+   ::
+
+      USE_SINGLE_PRECISION_JACOBIAN=TRUE
+
+   This can speed up the integration prevent the code from running out
+   of memory when run on GPUs.
+
+
+
 Integration errors
 ==================
 

diff --git a/integration/BackwardEuler/be_type.H b/integration/BackwardEuler/be_type.H
@@ -48,7 +48,7 @@ struct be_t {
     amrex::Real rtol_enuc;
 
     amrex::Array1D<amrex::Real, 1, int_neqs> y;
-    ArrayUtil::MathArray2D<1, int_neqs, 1, int_neqs> jac;
+    ArrayUtil::MathArray2D<jac_t, 1, int_neqs, 1, int_neqs> jac;
 
     short jacobian_type;
 };

diff --git a/integration/Make.package b/integration/Make.package
@@ -26,6 +26,9 @@ else
   CEXE_headers += integrator_setup_strang.H
 endif
 
+ifeq ($(USE_SINGLE_PRECISION_JACOBIAN), TRUE)
+  DEFINES += -DSINGLE_PRECISION_JACOBIAN
+endif
 
 ifeq ($(USE_ALL_NSE), TRUE)
   ifeq ($(USE_ALL_SDC), TRUE)

diff --git a/integration/VODE/vode_dvode.H b/integration/VODE/vode_dvode.H
@@ -170,14 +170,18 @@ int dvode (BurnT& state, DvodeT& vstate)
 
            if (vstate.n_step == 0) {
 #ifndef AMREX_USE_GPU
+           if (!state.suppress_failure_output) {
                std::cout << amrex::Font::Bold << amrex::FGColor::Red << "DVODE: too much accuracy requested at start of integration" << amrex::ResetDisplay << std::endl;
+           }
 #endif
                return IERR_TOO_MUCH_ACCURACY_REQUESTED;
            }
 
            // Too much accuracy requested for machine precision.
 #ifndef AMREX_USE_GPU
-           std::cout << "DVODE: too much accuracy requested" << std::endl;
+           if (!state.suppress_failure_output) {
+               std::cout << "DVODE: too much accuracy requested" << std::endl;
+           }
 #endif
            for (int i = 1; i <= int_neqs; ++i) {
                vstate.y(i) = vstate.yh(i,1);
@@ -197,7 +201,9 @@ int dvode (BurnT& state, DvodeT& vstate)
        if (kflag == -1) {
            // Error test failed repeatedly or with ABS(H) = HMIN.
 #ifndef AMREX_USE_GPU
-           std::cout << amrex::Font::Bold << amrex::FGColor::Red << "DVODE: error test failed repeatedly or with abs(H) = HMIN" << amrex::ResetDisplay << std::endl;
+           if (!state.suppress_failure_output) {
+               std::cout << amrex::Font::Bold << amrex::FGColor::Red << "DVODE: error test failed repeatedly or with abs(H) = HMIN" << amrex::ResetDisplay << std::endl;
+           }
 #endif
            // Set Y array, T, and optional output.
            for (int i = 1; i <= int_neqs; ++i) {
@@ -211,7 +217,9 @@ int dvode (BurnT& state, DvodeT& vstate)
        else if (kflag == -2) {
            // Convergence failed repeatedly or with ABS(H) = HMIN.
 #ifndef AMREX_USE_GPU
-           std::cout << amrex::Font::Bold << amrex::FGColor::Red << "DVODE: corrector convergence failed repeatedly or with abs(H) = HMIN" << amrex::ResetDisplay << std::endl;
+           if (!state.suppress_failure_output) {
+               std::cout << amrex::Font::Bold << amrex::FGColor::Red << "DVODE: corrector convergence failed repeatedly or with abs(H) = HMIN" << amrex::ResetDisplay << std::endl;
+           }
 #endif
            // Set Y array, T, and optional output.
            for (int i = 1; i <= int_neqs; ++i) {

diff --git a/integration/VODE/vode_type.H b/integration/VODE/vode_type.H
@@ -143,11 +143,11 @@ struct dvode_t
     amrex::Array1D<amrex::Real, 1, int_neqs> y;
 
     // Jacobian
-    ArrayUtil::MathArray2D<1, int_neqs, 1, int_neqs> jac;
+    ArrayUtil::MathArray2D<jac_t, 1, int_neqs, 1, int_neqs> jac;
 
 #ifdef ALLOW_JACOBIAN_CACHING
     // Saved Jacobian
-    ArrayUtil::MathArray2D<1, int_neqs, 1, int_neqs> jac_save;
+    ArrayUtil::MathArray2D<jac_t, 1, int_neqs, 1, int_neqs> jac_save;
 #endif
 
     // the Nordsieck history array

diff --git a/integration/integrator_data.H b/integration/integrator_data.H
@@ -52,6 +52,6 @@ constexpr int integrator_neqs ()
 
 using IArray1D = amrex::Array1D<short, 1, INT_NEQS>;
 using RArray1D = amrex::Array1D<amrex::Real, 1, INT_NEQS>;
-using RArray2D = ArrayUtil::MathArray2D<1, INT_NEQS, 1, INT_NEQS>;
+using RArray2D = ArrayUtil::MathArray2D<jac_t, 1, INT_NEQS, 1, INT_NEQS>;
 
 #endif
diff --git a/integration/integrator_setup_strang.H b/integration/integrator_setup_strang.H
@@ -181,7 +181,7 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state,
 
     // If we failed, print out the current state of the integration.
 
-    if (! state.success) {
+    if (! state.success && !state.suppress_failure_output) {
 #ifndef AMREX_USE_GPU
         std::cout << amrex::Font::Bold << amrex::FGColor::Red << "[ERROR] integration failed in net" << amrex::ResetDisplay << std::endl;
         std::cout << "istate = " << istate << std::endl;

diff --git a/integration/utils/circle_theorem.H b/integration/utils/circle_theorem.H
@@ -17,7 +17,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void circle_theorem_sprad(const amrex::Real time, BurnT& state, T& int_state, amrex::Real& sprad)
 {
 
-    ArrayUtil::MathArray2D<1, INT_NEQS, 1, INT_NEQS> jac_array;
+    ArrayUtil::MathArray2D<jac_t, 1, INT_NEQS, 1, INT_NEQS> jac_array;
 
     if (integrator_rp::jacobian == 1) {
         jac(time, state, int_state, jac_array);

diff --git a/interfaces/ArrayUtilities.H b/interfaces/ArrayUtilities.H
@@ -34,7 +34,7 @@ namespace ArrayUtil
         amrex::Real arr[(XHI-XLO+1)];
     };
 
-    template <int XLO, int XHI, int YLO, int YHI>
+    template <typename T, int XLO, int XHI, int YLO, int YHI>
     struct MathArray2D
     {
         AMREX_GPU_HOST_DEVICE AMREX_INLINE
@@ -71,7 +71,7 @@ namespace ArrayUtil
         }
 
         [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-        amrex::Real get (const int i, const int j) const noexcept {
+        T get (const int i, const int j) const noexcept {
             AMREX_ASSERT(i >= XLO && i <= XHI && j >= YLO && j <= YHI);
             return arr[i+j*(XHI-XLO+1)-(YLO*(XHI-XLO+1)+XLO)];
         }
@@ -84,18 +84,18 @@ namespace ArrayUtil
         }
 
         AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-        const amrex::Real& operator() (int i, int j) const noexcept {
+        const T& operator() (int i, int j) const noexcept {
             AMREX_ASSERT(i >= XLO && i <= XHI && j >= YLO && j <= YHI);
             return arr[i+j*(XHI-XLO+1)-(YLO*(XHI-XLO+1)+XLO)];
         }
 
         AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-        amrex::Real& operator() (int i, int j) noexcept {
+        T& operator() (int i, int j) noexcept {
             AMREX_ASSERT(i >= XLO && i <= XHI && j >= YLO && j <= YHI);
             return arr[i+j*(XHI-XLO+1)-(YLO*(XHI-XLO+1)+XLO)];
         }
 
-        amrex::Real arr[(XHI-XLO+1)*(YHI-YLO+1)];
+        T arr[(XHI-XLO+1)*(YHI-YLO+1)];
     };
 
     namespace Math

diff --git a/interfaces/burn_type.H b/interfaces/burn_type.H
@@ -56,7 +56,15 @@ const int SVAR  = SEDEN+1;
 const int SVAR_EVOLVE = SFX;
 
 // this is the data type of the dense Jacobian that the network wants.
-using  JacNetArray2D = ArrayUtil::MathArray2D<1, neqs, 1, neqs>;
+
+// the Jacobian can have a different type that our system
+#ifdef SINGLE_PRECISION_JACOBIAN
+using jac_t = float;
+#else
+using jac_t = amrex::Real;
+#endif
+
+using  JacNetArray2D = ArrayUtil::MathArray2D<jac_t, 1, neqs, 1, neqs>;
 
 struct burn_t
 {
@@ -154,6 +162,9 @@ struct burn_t
   // integrator error code
   short error_code{};
 
+  // suppress verbose failure diagnostics (useful for GPU retries)
+  bool suppress_failure_output{};
+
 };
 
 

diff --git a/networks/rhs.H b/networks/rhs.H
@@ -1357,7 +1357,7 @@ void rhs (burn_t& burn_state, amrex::Array1D<amrex::Real, 1, nrhs>& ydot)
 
 // Analytical Jacobian
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void jac (burn_t& burn_state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac)
+void jac (burn_t& burn_state, ArrayUtil::MathArray2D<jac_t, 1, neqs, 1, neqs>& jac)
 {
     rhs_state_t<autodiff::dual> rhs_state;
 
@@ -1521,7 +1521,7 @@ void actual_rhs (burn_t& state, amrex::Array1D<amrex::Real, 1, neqs>& ydot)
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void actual_jac (burn_t& state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac)
+void actual_jac (burn_t& state, ArrayUtil::MathArray2D<jac_t, 1, neqs, 1, neqs>& jac)
 {
     RHS::jac(state, jac);
 }

diff --git a/unit_test/burn_cell_primordial_chem/GNUmakefile b/unit_test/burn_cell_primordial_chem/GNUmakefile
@@ -2,7 +2,7 @@ PRECISION  = DOUBLE
 PROFILE    = FALSE
 
 # Set DEBUG to TRUE if debugging
-DEBUG      = TRUE
+DEBUG      = FALSE
 
 DIM        = 1
 
@@ -15,7 +15,7 @@ USE_CUDA   = FALSE
 USE_REACT = TRUE
 
 # Set USE_MICROPHYSICS_DEBUG to TRUE if debugging
-USE_MICROPHYSICS_DEBUG = TRUE
+USE_MICROPHYSICS_DEBUG = FALSE
 
 EBASE = main