[WIP] allow for an alternate EOS table for the helmholtz EOS

.github/workflows/good_defines.txt

@@ -27,6 +27,7 @@ SIMPLIFIED_SDC
 SINGLE_PRECISION_JACOBIAN
 STRANG
 TRUE_SDC
+USE_NEW_HELMHOLTZ_TABLE
 _OPENMP
 _WIN32
 __cplusplus

EOS/helmholtz/Make.package

@@ -1,3 +1,4 @@
 CEXE_headers += actual_eos.H
 CEXE_headers += actual_eos_data.H
 CEXE_sources += actual_eos_data.cpp
+CEXE_headers += mp_math.H

EOS/helmholtz/actual_eos.H

@@ -19,6 +19,9 @@
 #include <actual_eos_data.H>
 #include <cmath>
 #include <vector>
+#if defined(USE_NEW_HELMHOLTZ_TABLE)
+#include <mp_math.H>
+#endif
 
 // Frank Timmes Helmholtz based Equation of State
 // https://cococubed.com/
@@ -161,9 +164,15 @@ void apply_electrons (T& state)
     amrex::Real din = state.y_e * state.rho;
 
     // hash locate this temperature and density
+#if defined(USE_NEW_HELMHOLTZ_TABLE)
+    // we use the fastmath variant of log
+    int jat = int((mp::fastlg2(state.T) - tlo) * tstpi) + 1;
+    int iat = int((mp::fastlg2(din) - dlo) * dstpi) + 1;
+#else
     int jat = int((std::log10(state.T) - tlo) * tstpi) + 1;
-    jat = amrex::Clamp(jat, 1, jmax-1) - 1;
     int iat = int((std::log10(din) - dlo) * dstpi) + 1;
+#endif
+    jat = amrex::Clamp(jat, 1, jmax-1) - 1;
     iat = amrex::Clamp(iat, 1, imax-1) - 1;
 
     amrex::Real fi[36];
@@ -1312,10 +1321,19 @@ void actual_eos_init ()
 
     for (int j = 0; j < jmax; ++j) {
         amrex::Real tsav = tlo + j * tstp;
+#if defined(USE_NEW_HELMHOLTZ_TABLE)
+        t[j] = mp::fastpow2(tsav);
+#else
         t[j] = std::pow(10.0e0_rt, tsav);
+#endif
+
         for (int i = 0; i < imax; ++i) {
             amrex::Real dsav = dlo + i * dstp;
+#if defined(USE_NEW_HELMHOLTZ_TABLE)
+            d[i] = mp::fastpow2(dsav);
+#else
             d[i] = std::pow(10.0e0_rt, dsav);
+#endif
         }
     }
 
@@ -1333,7 +1351,12 @@ void actual_eos_init ()
 
         // open the table
         std::ifstream table;
+
+#if defined(USE_NEW_HELMHOLTZ_TABLE)
+        table.open("helm_table_fastmath_p128_q800.dat");
+#else
         table.open("helm_table.dat");
+#endif
 
         if (!table.is_open()) {
             // the table was not present or we could not open it; abort
@@ -1462,10 +1485,18 @@ void actual_eos_init ()
 
     // Set up the minimum and maximum possible densities.
 
+#if defined(USE_NEW_HELMHOLTZ_TABLE)
+    EOSData::mintemp = mp::fastpow2(tlo);
+    EOSData::maxtemp = mp::fastpow2(thi);
+    EOSData::mindens = mp::fastpow2(dlo);
+    EOSData::maxdens = mp::fastpow2(dhi);
+#else
     EOSData::mintemp = std::pow(10.e0_rt, tlo);
     EOSData::maxtemp = std::pow(10.e0_rt, thi);
     EOSData::mindens = std::pow(10.e0_rt, dlo);
     EOSData::maxdens = std::pow(10.e0_rt, dhi);
+#endif
+
 }
 
 

EOS/helmholtz/actual_eos_data.H

@@ -4,6 +4,10 @@
 #include <AMReX.H>
 #include <AMReX_REAL.H>
 
+#if defined(USE_NEW_HELMHOLTZ_TABLE)
+#include <fundamental_constants.H>
+#endif
+
 namespace helmholtz
 {
 
@@ -14,6 +18,22 @@ namespace helmholtz
 
     // for the tables
 
+#if defined(USE_NEW_HELMHOLTZ_TABLE)
+    constexpr int imax = 421;
+    constexpr int jmax = 161;
+
+    // these are log_2 using the fastlog approximation
+    constexpr amrex::Real tlo   = 9.953125_rt;
+    constexpr amrex::Real thi   = 36.4551915228366851806640625_rt;
+    constexpr amrex::Real tstp  = (thi - tlo) / (static_cast<amrex::Real>(jmax-1));
+    constexpr amrex::Real tstpi = 1.0e0_rt / tstp;
+    constexpr amrex::Real dlo   = -33.2820130816_rt;
+    constexpr amrex::Real dhi   = 36.4551915228366851806640625_rt;
+    constexpr amrex::Real dstp  = (dhi - dlo) / (static_cast<amrex::Real>(imax-1));
+    constexpr amrex::Real dstpi = 1.0e0_rt / dstp;
+
+#else
+
     constexpr int imax = 541;
     constexpr int jmax = 201;
 
@@ -25,6 +45,7 @@ namespace helmholtz
     constexpr amrex::Real dhi   = 15.0e0_rt;
     constexpr amrex::Real dstp  = (dhi - dlo) / (static_cast<amrex::Real>(imax-1));
     constexpr amrex::Real dstpi = 1.0e0_rt / dstp;
+#endif
 
     extern AMREX_GPU_MANAGED amrex::Real d[imax];
     extern AMREX_GPU_MANAGED amrex::Real t[jmax];
@@ -55,11 +76,19 @@ namespace helmholtz
     extern AMREX_GPU_MANAGED amrex::Real ddi_sav[imax];
     extern AMREX_GPU_MANAGED amrex::Real dd2i_sav[imax];
 
+#if defined(USE_NEW_HELMHOLTZ_TABLE)
+    constexpr amrex::Real h = C::hplanck;
+    constexpr amrex::Real avo_eos = C::n_A;
+    constexpr amrex::Real kerg = C::k_B;
+    constexpr amrex::Real amu = C::m_u;
+#else
     // 2006 CODATA physical constants
     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
+
 }
 
 #endif