generalize in_nse() to work with an eos_t (#1424)

AMReX-Astro · Dec 29, 2023 · fd59277 · fd59277
1 parent 0c3a1d0
commit fd59277
Showing 1 changed file with 47 additions and 19 deletions.
diff --git a/nse_tabular/nse_table_check.H b/nse_tabular/nse_table_check.H
@@ -14,17 +14,18 @@
 #include <extern_parameters.H>
 
 #include <burn_type.H>
-#include <eos.H>
+#include <eos_type.H>
 
 
+template <typename T>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-bool in_nse(burn_t& state, const bool relax = false) {
+bool in_nse(T& state, const bool relax = false) {
 
     using namespace Species;
 
     Real f = relax ? nse_relax_factor : 1.0_rt;
 
-    state.nse = false;
+    bool nse_check = false;
 
     if (state.rho > f * rho_nse && state.T > f * T_nse) {
 
@@ -33,33 +34,60 @@ bool in_nse(burn_t& state, const bool relax = false) {
         // and He-group (for us, that is H1, He3, He4)
         // also make sure there is not a lot of O16 (O16 burning will dominate then) or Si28
 
+        Real Fe_group{};
+        Real C_group{};
+        Real He_group{};
+        Real O_group{};
+        Real Si_group{};
+
+        if constexpr (std::is_same<T, burn_t>::value) {
+            // for a burn_t, we need to use a different field
+            // depending on whether we are strang of simplified-SDC
+
 #ifdef STRANG
-        Real Fe_group = state.xn[Cr48-1] + state.xn[Fe52-1] +
-                        state.xn[Fe54-1] + state.xn[Ni56-1];
-        Real C_group = state.xn[C12-1] + state.xn[N14-1];
-        Real He_group = state.xn[H1-1] + state.xn[He3-1] + state.xn[He4-1];
-        Real O_group = state.xn[O16-1];
-        Real Si_group = state.xn[Si28-1];
+            Fe_group = state.xn[Cr48-1] + state.xn[Fe52-1] +
+                       state.xn[Fe54-1] + state.xn[Ni56-1];
+            C_group = state.xn[C12-1] + state.xn[N14-1];
+            He_group = state.xn[H1-1] + state.xn[He3-1] + state.xn[He4-1];
+            O_group = state.xn[O16-1];
+            Si_group = state.xn[Si28-1];
 #else
-        // we need to get the mass fractions from the conserved state passed in
-
-        Real Fe_group = (state.y[SFS+Cr48-1] + state.y[SFS+Fe52-1] +
-                         state.y[SFS+Fe54-1] + state.y[SFS+Ni56-1]) / state.rho;
-        Real C_group = (state.y[SFS+C12-1] + state.y[SFS+N14-1]) / state.rho;
-        Real He_group = (state.y[SFS+H1-1] + state.y[SFS+He3-1] + state.y[SFS+He4-1]) / state.rho;
-        Real O_group = state.y[SFS+O16-1] / state.rho;
-        Real Si_group = state.y[SFS+Si28-1] / state.rho;
+            // we need to get the mass fractions from the conserved state passed in
+
+            Fe_group = (state.y[SFS+Cr48-1] + state.y[SFS+Fe52-1] +
+                        state.y[SFS+Fe54-1] + state.y[SFS+Ni56-1]) / state.rho;
+            C_group = (state.y[SFS+C12-1] + state.y[SFS+N14-1]) / state.rho;
+            He_group = (state.y[SFS+H1-1] + state.y[SFS+He3-1] + state.y[SFS+He4-1]) / state.rho;
+            O_group = state.y[SFS+O16-1] / state.rho;
+            Si_group = state.y[SFS+Si28-1] / state.rho;
 #endif
 
+        } else {
+
+            // this covers the various eos_t's
+            Fe_group = state.xn[Cr48-1] + state.xn[Fe52-1] +
+                       state.xn[Fe54-1] + state.xn[Ni56-1];
+            C_group = state.xn[C12-1] + state.xn[N14-1];
+            He_group = state.xn[H1-1] + state.xn[He3-1] + state.xn[He4-1];
+            O_group = state.xn[O16-1];
+            Si_group = state.xn[Si28-1];
+        }
+
         if (Fe_group + He_group > f * He_Fe_nse &&
             C_group < C_nse / f &&
             O_group < O_nse / f &&
             Si_group < Si_nse / f) {
-            state.nse = true;
+
+            nse_check = true;
         }
+
+    }
+
+    if constexpr (std::is_same<T, burn_t>::value) {
+        state.nse = nse_check;
     }
 
-    return state.nse;
+    return nse_check;
 
 }