Merge branch 'development' into xrb_iginition

Author: zhichen3

.github/workflows/nse_net.yml

@@ -65,3 +65,23 @@ jobs:
         run: |
           cd unit_test/nse_net_cell/make_table
           diff -I "^Initializing AMReX" -I "^AMReX" -I "^reading in reaclib rates" test.out ci-benchmarks/nse_net_make_table_unit_test.out
+
+      - name: Compile, nse_compatibility (NSE_NET, ase-test, integration)
+        run: |
+          cd nse_solver/nse_compatibility
+          make realclean
+          make NETWORK_DIR=he-burn/ase-test -j 4
+
+      - name: Run nse_compatibility (NSE_NET, ase-test, integration)
+        run: |
+          cd nse_solver/nse_compatibility
+          ./main3d.gnu.ex inputs_nse_compatibility amrex.fpe_trap_{invalid,zero,overflow}=1 > test.out
+
+      - name: Print backtrace
+        if: ${{ failure() && hashFiles('nse_solver/nse_compatibility/Backtrace.0') != '' }}
+        run: cat nse_solver/nse_compatibility/Backtrace.0
+
+      - name: Compare to stored output (NSE_NET, ase-test, integration)
+        run: |
+          cd nse_solver/nse_compatibility
+          diff -I "^Initializing AMReX" -I "^AMReX" -I "^reading in reaclib rates" test.out ci-benchmarks/nse_net_integration_unit_test.out

Docs/source/ode_integrators.rst

@@ -248,7 +248,7 @@ fractions.  Looser than this can produce large errors.
 Controlling Species $\sum_k X_k = 1$
 ====================================
 
-.. index:: integrator.renormalize_abundances, integrator.SMALL_X_SAFE, integrator.do_species_clip
+.. index:: integrator.renormalize_abundances, integrator.SMALL_X_SAFE, integrator.do_species_clip, integrator.do_corrector_validation
 
 The ODE integrators don't know about the constraint that
 
@@ -275,10 +275,18 @@ constraint on the intermediate states during the integration.
   The default is ``1.e-30``.
 
 * ``integrator.do_species_clip`` : this enforces that the mass fractions
-  all in $[\mathtt{SMALL\_X\_SAFE}, 1.0]$.
-
-  This is enabled by default.
-
+  all in $[\mathtt{SMALL\_X\_SAFE}, 1.0]$ before calling the network righthand
+  side function.
+
+  This is off by default.  Turning this on can sometimes make the integrator
+  work a lot harder.
+
+* ``integrator.do_corrector_validation`` : in the nonlinear solve
+  corrector loop, when we get a corrected integration state, do we
+  check to make sure the mass fractions are valid before calling the
+  righthand function?  This is needed in some cases if
+  ``integrator.do_species_clip`` is disabled.  Note: this is not
+  implemented for every integrator.
 
 
 Retry Mechanism

EOS/helmholtz/actual_eos.H

@@ -6,6 +6,8 @@
 #include <iostream>
 #include <fstream>
 #include <sstream>
+#include <numbers>
+
 #include <AMReX.H>
 #include <AMReX_REAL.H>
 #include <AMReX_ParallelDescriptor.H>

integration/VODE/vode_dvnlsd.H

@@ -145,6 +145,36 @@ amrex::Real dvnlsd (int& NFLAG, BurnT& state, DvodeT& vstate)
                 vstate.y(i) = vstate.yh(i,1) + vstate.acor(i);
             }
 
+            // sometime VODE goes way off tangent.  If these mass fractions
+            // are really bad, then let's just bail now
+
+            if (integrator_rp::do_corrector_validation && !integrator_rp::use_number_densities) {
+#ifdef SDC
+                const amrex::Real rho_current = state.rho_orig + vstate.tn * state.ydot_a[SRHO];
+                const amrex::Real thresh = species_failure_tolerance * rho_current;
+#else
+                const amrex::Real thresh = species_failure_tolerance;
+#endif
+                bool fail{};
+                for (int i = 1; i <= NumSpec; ++i) {
+                    if (vstate.y(i) < -thresh) {
+                        fail = true;
+                        break;
+                    }
+                }
+
+                // this resets the flags in the same fashion as is done above
+                // for a singular matrix.  ICF = 2 means "unrecoverable error",
+                // IPUP = 1 forces a Jacobian reevaluation.  The return value,
+                // ACNRM is ignored in this case (since NFLAG != 0)
+                if (fail) {
+                    NFLAG = -1;
+                    vstate.ICF = 2;
+                    vstate.IPUP = 1;
+                    return ACNRM;
+                }
+            }
+
             // Test for convergence. If M > 0, an estimate of the convergence
             // rate constant is stored in CRATE, and this is used in the test.
 

integration/_parameters

@@ -91,7 +91,11 @@ retry_atol_enuc                real      -1
 
 # in the clean_state process, do we clip the species such that they
 # are in [0, 1]?
-do_species_clip              bool         1
+do_species_clip              bool         0
+
+# in the corrector loop, do we check if the predicted state is
+# valid (X > 0) before calling the RHS?
+do_corrector_validation      bool         1
 
 # flag for turning on the use of number densities for all species
 use_number_densities     bool     0

networks/he-burn/ase-test/Make.package

@@ -0,0 +1,17 @@
+CEXE_headers += network_properties.H
+
+ifeq ($(USE_REACT),TRUE)
+  CEXE_sources += actual_network_data.cpp
+  CEXE_headers += actual_network.H
+  CEXE_headers += tfactors.H
+  CEXE_headers += interp_tools.H
+  CEXE_headers += partition_functions.H
+  CEXE_sources += partition_functions_data.cpp
+  CEXE_headers += actual_rhs.H
+  CEXE_headers += reaclib_rates.H
+  CEXE_headers += table_rates.H
+  CEXE_sources += table_rates_data.cpp
+  CEXE_headers += temperature_table_rates.H
+  USE_SCREENING = TRUE
+  USE_NEUTRINOS = TRUE
+endif

networks/he-burn/ase-test/_parameters

@@ -0,0 +1,4 @@
+@namespace: network
+
+disable_p_C12_to_N13_reaclib    int     0
+disable_He4_N13_to_p_O16_reaclib    int     0