Merge branch 'nse_neutrinos' into nse_second_order_refactor

AMReX-Astro · Nov 22, 2023 · 6fa7ac5 · 6fa7ac5
2 parents 58b85f1 + d42279e
commit 6fa7ac5
Show file tree

Hide file tree

Showing 19 changed files with 8,253 additions and 42 deletions.
diff --git a/EOS/helmholtz/actual_eos.H b/EOS/helmholtz/actual_eos.H
@@ -155,9 +155,9 @@ void apply_electrons (T& state)
 
     // hash locate this temperature and density
     int jat = int((std::log10(state.T) - tlo) * tstpi) + 1;
-    jat = amrex::max(1, amrex::min(jat, jtmax-1)) - 1;
+    jat = amrex::max(1, amrex::min(jat, jmax-1)) - 1;
     int iat = int((std::log10(din) - dlo) * dstpi) + 1;
-    iat = amrex::max(1, amrex::min(iat, itmax-1)) - 1;
+    iat = amrex::max(1, amrex::min(iat, imax-1)) - 1;
 
     Real fi[36];
 
@@ -477,9 +477,9 @@ void apply_ions (T& state)
 {
     using namespace helmholtz;
 
-    const Real pi      = 3.1415926535897932384e0_rt;
-    const Real sioncon = (2.0e0_rt * pi * amu * kerg)/(h*h);
-    const Real kergavo = kerg * avo_eos;
+    constexpr Real pi      = 3.1415926535897932384e0_rt;
+    constexpr Real sioncon = (2.0e0_rt * pi * amu * kerg)/(h*h);
+    constexpr Real kergavo = kerg * avo_eos;
 
     Real deni = 1.0e0_rt / state.rho;
     Real tempi = 1.0e0_rt / state.T;
@@ -553,14 +553,14 @@ void apply_radiation (T& state)
 {
     using namespace helmholtz;
 
-    const Real clight  = 2.99792458e10_rt;
+    constexpr Real clight  = 2.99792458e10_rt;
 #ifdef RADIATION
-    const Real ssol    = 0.0e0_rt;
+    constexpr Real ssol    = 0.0e0_rt;
 #else
-    const Real ssol    = 5.67051e-5_rt;
+    constexpr Real ssol    = 5.67051e-5_rt;
 #endif
-    const Real asol    = 4.0e0_rt * ssol / clight;
-    const Real asoli3  = asol/3.0e0_rt;
+    constexpr Real asol    = 4.0e0_rt * ssol / clight;
+    constexpr Real asoli3  = asol/3.0e0_rt;
 
     Real deni = 1.0e0_rt / state.rho;
     Real tempi = 1.0e0_rt / state.T;
@@ -638,19 +638,19 @@ void apply_coulomb_corrections (T& state)
     using namespace helmholtz;
 
     // Constants used for the Coulomb corrections
-    const Real a1 = -0.898004e0_rt;
-    const Real b1 =  0.96786e0_rt;
-    const Real c1 =  0.220703e0_rt;
-    const Real d1 = -0.86097e0_rt;
-    const Real e1 =  2.5269e0_rt;
-    const Real a2 =  0.29561e0_rt;
-    const Real b2 =  1.9885e0_rt;
-    const Real c2 =  0.288675e0_rt;
-    const Real qe   = 4.8032042712e-10_rt;
-    const Real esqu = qe * qe;
-    const Real onethird = 1.0e0_rt/3.0e0_rt;
-    const Real forth = 4.0e0_rt/3.0e0_rt;
-    const Real pi    = 3.1415926535897932384e0_rt;
+    constexpr Real a1 = -0.898004e0_rt;
+    constexpr Real b1 =  0.96786e0_rt;
+    constexpr Real c1 =  0.220703e0_rt;
+    constexpr Real d1 = -0.86097e0_rt;
+    constexpr Real e1 =  2.5269e0_rt;
+    constexpr Real a2 =  0.29561e0_rt;
+    constexpr Real b2 =  1.9885e0_rt;
+    constexpr Real c2 =  0.288675e0_rt;
+    constexpr Real qe   = 4.8032042712e-10_rt;
+    constexpr Real esqu = qe * qe;
+    constexpr Real onethird = 1.0e0_rt/3.0e0_rt;
+    constexpr Real forth = 4.0e0_rt/3.0e0_rt;
+    constexpr Real pi    = 3.1415926535897932384e0_rt;
 
     [[maybe_unused]] Real pcoul    = 0.e0_rt;
     [[maybe_unused]] Real dpcouldd = 0.e0_rt;
@@ -1242,7 +1242,7 @@ void actual_eos (I input, T& state)
 
     using namespace helmholtz;
 
-    const int max_newton = 100;
+    constexpr int max_newton = 100;
 
     bool single_iter, converged;
     int var{}, dvar{}, var1, var2;
@@ -1306,7 +1306,6 @@ void actual_eos_init ()
 
     Real dth, dt2, dti, dt2i;
     Real dd, dd2, ddi, dd2i;
-    Real tsav, dsav;
 
     // Read in the runtime parameters
 
@@ -1316,8 +1315,6 @@ void actual_eos_init ()
     dtol = eos_dtol;
 
     //    read the helmholtz free energy table
-    itmax = imax;
-    jtmax = jmax;
     tlo   = 3.0e0_rt;
     thi   = 13.0e0_rt;
     tstp  = (thi - tlo) / ((Real) (jmax-1));
@@ -1328,10 +1325,10 @@ void actual_eos_init ()
     dstpi = 1.0e0_rt / dstp;
 
     for (int j = 0; j < jmax; ++j) {
-        tsav = tlo + j * tstp;
+        Real tsav = tlo + j * tstp;
         t[j] = std::pow(10.0e0_rt, tsav);
         for (int i = 0; i < imax; ++i) {
-            dsav = dlo + i * dstp;
+            Real dsav = dlo + i * dstp;
             d[i] = std::pow(10.0e0_rt, dsav);
         }
     }
@@ -1508,4 +1505,3 @@ bool is_input_valid (I input)
 }
 
 #endif
-
diff --git a/EOS/helmholtz/actual_eos_data.H b/EOS/helmholtz/actual_eos_data.H
@@ -6,13 +6,13 @@
 
 namespace helmholtz
 {
-    extern AMREX_GPU_MANAGED bool do_coulomb;
-    extern AMREX_GPU_MANAGED bool input_is_constant;
+    extern AMREX_GPU_MANAGED int do_coulomb;
+    extern AMREX_GPU_MANAGED int input_is_constant;
 
     // for the tables
 
-    const int imax = 541;
-    const int jmax = 201;
+    constexpr int imax = 541;
+    constexpr int jmax = 201;
 
     extern AMREX_GPU_MANAGED int itmax;
     extern AMREX_GPU_MANAGED int jtmax;
@@ -57,10 +57,10 @@ namespace helmholtz
     extern AMREX_GPU_MANAGED amrex::Real dd2i_sav[imax];
 
     // 2006 CODATA physical constants
-    const amrex::Real h = 6.6260689633e-27;
-    const amrex::Real avo_eos = 6.0221417930e23;
-    const amrex::Real kerg = 1.380650424e-16;
-    const amrex::Real amu = 1.66053878283e-24;
+    constexpr amrex::Real h = 6.6260689633e-27;
+    constexpr amrex::Real avo_eos = 6.0221417930e23;
+    constexpr amrex::Real kerg = 1.380650424e-16;
+    constexpr amrex::Real amu = 1.66053878283e-24;
 }
 
 #endif
diff --git a/EOS/helmholtz/actual_eos_data.cpp b/EOS/helmholtz/actual_eos_data.cpp
@@ -1,7 +1,7 @@
 #include <actual_eos_data.H>
 
-AMREX_GPU_MANAGED bool helmholtz::do_coulomb;
-AMREX_GPU_MANAGED bool helmholtz::input_is_constant;
+AMREX_GPU_MANAGED int helmholtz::do_coulomb;
+AMREX_GPU_MANAGED int helmholtz::input_is_constant;
 
 AMREX_GPU_MANAGED int helmholtz::itmax;
 AMREX_GPU_MANAGED int helmholtz::jtmax;

diff --git a/networks/subch_base/Make.package b/networks/subch_base/Make.package
@@ -0,0 +1,14 @@
+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 += partition_functions.H
+  CEXE_headers += actual_rhs.H
+  CEXE_headers += reaclib_rates.H
+  CEXE_headers += table_rates.H
+  CEXE_sources += table_rates_data.cpp
+  USE_SCREENING = TRUE
+  USE_NEUTRINOS = TRUE
+endif
diff --git a/networks/subch_base/README.md b/networks/subch_base/README.md
@@ -0,0 +1,8 @@
+# subch_base
+
+This network is a further simplification of subch_simple by
+excluding N14, F18, Ne21, and Na22.
+These nuclei were added to include n14(a, g)f18(a, p)ne21 rate sequence,
+which allows an enhancement to the c12(p, g)n13(a, p)o16 rate due to the
+proton release. However, we find that the effect is not very significant,
+and they also caused an artificial endpoint at Na22.
diff --git a/networks/subch_base/_parameters b/networks/subch_base/_parameters
@@ -0,0 +1,4 @@
+@namespace: network
+
+disable_p_C12_to_N13    int     0
+disable_He4_N13_to_p_O16    int     0