Merge pull request #20 from lanl/dholladay00/optimized_logs

Adding optimized log implementation using frexp and approximation to …

.gitmodules

@@ -10,9 +10,6 @@
 [submodule "utils/variant"]
 	path = utils/variant
 	url = https://github.com/mpark/variant.git
-[submodule "utils/herumi-fmath"]
-	path = utils/herumi-fmath
-	url = https://github.com/herumi/fmath.git
 [submodule "utils/kokkos"]
 	path = utils/kokkos
 	url = https://github.com/kokkos/kokkos.git

CMakeLists.txt

@@ -35,6 +35,10 @@ option (SINGULARITY_BUILD_CLOSURE "Mixed cell closure" ON)
 option (SINGULARITY_TEST_SESAME "Test the Sesame table readers" OFF)
 option (SINGULARITY_TEST_STELLAR_COLLAPSE "Test the stellar collapse table readers" OFF)
 option (SINGULARITY_EOS_SKIP_EXTRAP "Turn off eospac extrap checks" OFF)
+option (SINGULARITY_USE_SINGLE_LOGS "Use single precision logs. Can harm accuracy." OFF)
+option (SINGULARITY_FMATH_USE_ORDER_4 "4th order interpolant for fast logs. Default is 7th order." OFF)
+option (SINGULARITY_FMATH_USE_ORDER_5 "5th order interpolant for fast logs. Default is 7th order." OFF)
+
 if (SINGULARITY_SUBMODULE_MODE)
   set(SINGULARITY_BETTER_DEBUG_FLAGS OFF CACHE BOOL "" FORCE)
   set(SINGULARITY_HIDE_MORE_WARNINGS ON CACHE BOOL "" FORCE)
@@ -280,6 +284,23 @@ if (SINGULARITY_USE_EOSPAC)
   target_link_libraries(singularity-eos::libs INTERFACE EOSPAC::eospac)
 endif ()
 
+# fast math
+if (SINGULARITY_USE_SINGLE_LOGS)
+  target_compile_definitions(singularity-eos::flags INTERFACE
+      SINGULARITY_USE_SINGLE_LOGS)
+endif()
+if (SINGULARITY_FMATH_USE_ORDER_4 AND SINGULARITY_FMATH_USE_ORDER_5)
+  message(ERROR "Order 4 and order 5 interpolation both specified. Please specify only one or zero. If no specification is made, order 7 will be used.")
+endif()
+if (SINGULARITY_FMATH_USE_ORDER_4)
+  target_compile_definitions(singularity-eos::flags INTERFACE
+    SINGULAIRTY_FMATH_USE_ORDER_4)
+endif()
+if (SINGULARITY_FMATH_USE_ORDER_5)
+  target_compile_definitions(singularity-eos::flags INTERFACE
+    SINGULAIRTY_FMATH_USE_ORDER_5)
+endif()
+
 # required package ports of call
 find_package(PortsofCall REQUIRED)
 target_link_libraries(singularity-eos::flags INTERFACE PortsofCall::PortsofCall)
@@ -348,10 +369,6 @@ install(DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/utils/fast-math"
   DESTINATION "include"
   FILES_MATCHING PATTERN "*.hpp"
   )
-install(DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/utils/herumi-fmath"
-  DESTINATION "include"
-  FILES_MATCHING PATTERN "*.hpp"
-  )
 install(DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/utils/root-finding-1d"
   DESTINATION "include"
   FILES_MATCHING PATTERN "*.hpp"

sesame2spiner/examples/airtest.json

@@ -0,0 +1,32 @@
+//======================================================================
+// Example input deck for sesame2spiner,
+// a tool for converting eospac to spiner
+// Author: Jonah Miller ([email protected])
+// © 2021. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S.  Department of Energy/National
+// Nuclear Security Administration. All rights in the program are
+// reserved by Triad National Security, LLC, and the U.S. Department of
+// Energy/National Nuclear Security Administration. The Government is
+// granted for itself and others acting on its behalf a nonexclusive,
+// paid-up, irrevocable worldwide license in this material to reproduce,
+// prepare derivative works, distribute copies to the public, perform
+// publicly and display publicly, and to permit others to do so.
+//======================================================================
+{
+    "savename"  : "materials.sp5",
+    "materials" : [
+  {
+      "matid" : 5030,
+      "name"  : "air",
+      /* These set the number of grid points per decade
+         for each variable. The default is 50 points
+         per decade.
+      */
+      "numrho/decade" : 40,
+      "numT/decade"   : 40,
+      "numSie/decade" : 40
+  }
+    ]
+}

sesame2spiner/examples/materials.json

@@ -34,7 +34,7 @@
       "matid" : 2961,
       "name"  : "titanium",
       /* These set the number of grid points per decade
-         for each variable. The default is 20 points
+         for each variable. The default is 50 points
          per decade.
       */
       "numrho/decade" : 30,

sesame2spiner/generate_files.cpp

@@ -302,6 +302,15 @@ void getMatBounds(int i,
   Real rhoAnchor = metadata.normalDensity;
   Real TAnchor = 298.15;
 
+  // Forces density and temperature to be in a region where an offset
+  // is not needed. This improves resolution at low densities and
+  // temperatures.
+
+  // Extrapolation and other resolution tricks will be explored in the
+  // future.
+  if (rhoMin < STRICTLY_POS_MIN) rhoMin = STRICTLY_POS_MIN;
+  if (TMin < STRICTLY_POS_MIN) TMin = STRICTLY_POS_MIN;
+
   lRhoBounds = Bounds(rhoMin,rhoMax,numRho,true,shrinklRhoBounds,rhoAnchor);
   lTBounds   = Bounds(TMin,TMax,numT,true,shrinklTBounds,TAnchor);
   leBounds   = Bounds(sieMin,sieMax,numSie,true,shrinkleBounds);

sesame2spiner/generate_files.hpp

@@ -27,6 +27,7 @@
 
 using nlohmann::json;
 constexpr int PPD_DEFAULT = 50;
+constexpr Real STRICTLY_POS_MIN = 1e-9;
 
 herr_t saveMaterial(hid_t loc,
                     const SesameMetadata& metadata,

sesame2spiner/io_eospac.hpp

@@ -74,28 +74,23 @@ class Bounds {
     if (convertToLog) {
       // should be single-precision epsilon b/c that's whats used for the logs
       constexpr Real epsilon = std::numeric_limits<float>::epsilon();
-      const Real min_offset = std::max(10*std::abs(epsilon),
-                                       std::abs(epsilon*max));
+      const Real min_offset = 10*std::abs(epsilon);
       if (min < 0) offset = std::abs(min) + min_offset;
       else if ( min == 0 ) {
         offset = min_offset;
       }
       min += offset;
       max += offset;
-      // min = std::max(std::max(10*std::abs(Spiner::EPS),
-      //                         std::abs(Spiner::EPS*max)),
-      //                min);
-      // if (min <= 0) min = std::abs(10*Spiner::EPS);
-      // if (min <= 0) min = std::abs(Spiner::EPS);
-      min = BDMath::log10(std::abs(min)); // fast, floating-point log
-      max = BDMath::log10(std::abs(max));
+
+      min = std::log10(std::abs(min));
+      max = std::log10(std::abs(max));
       Real delta = max - min;
       min += 0.5*shrinkRange*delta;
       max -= 0.5*shrinkRange*delta;
 
       if (!(std::isnan(anchor_point))) {
         anchor_point += offset;
-        anchor_point = BDMath::log10(std::abs(anchor_point));
+        anchor_point = std::log10(std::abs(anchor_point));
       }
     }
 

sesame2spiner/parser.hpp

@@ -42,7 +42,7 @@ const std::string EXAMPLESTRING = R"(
       "matid" : 2961,
       "name"  : "titanium",
       /* These set the number of grid points per decade
-         for each variable. The default is 20 points
+         for each variable. The default is 50 points
          per decade.
       */
       "numrho/decade" : 30,

singularity-eos/eos/eos.hpp

@@ -869,7 +869,7 @@ class SpinerEOSDependsRhoT {
   toLog_(const Real x, const Real offset) const noexcept {
     // return std::log10(x + offset + EPS);
     // return std::log10(std::abs(std::max(x,-offset) + offset)+EPS);
-    return BDMath::log10(std::abs(std::max(x, -offset) + offset) + EPS);
+    return Math::log10(std::abs(std::max(x, -offset) + offset) + EPS);
   }
   PORTABLE_INLINE_FUNCTION Real __attribute__((always_inline))
   fromLog_(const Real lx, const Real offset) const noexcept {
@@ -1076,7 +1076,7 @@ class SpinerEOSDependsRhoSie {
 
   PORTABLE_INLINE_FUNCTION Real toLog_(const Real x, const Real offset) const {
     // return std::log10(std::abs(std::max(x,-offset) + offset)+EPS);
-    return BDMath::log10(std::abs(std::max(x, -offset) + offset) + EPS);
+    return Math::log10(std::abs(std::max(x, -offset) + offset) + EPS);
   }
   PORTABLE_INLINE_FUNCTION Real fromLog_(const Real lx,
                                          const Real offset) const {
@@ -1262,7 +1262,7 @@ class StellarCollapse {
   toLog_(const Real x, const Real offset) const noexcept {
     // return std::log10(x + offset + EPS);
     // return std::log10(std::abs(std::max(x,-offset) + offset)+EPS);
-    return BDMath::log10(std::abs(std::max(x, -offset) + offset) + EPS);
+    return Math::log10(std::abs(std::max(x, -offset) + offset) + EPS);
   }
   PORTABLE_INLINE_FUNCTION Real __attribute__((always_inline))
   fromLog_(const Real lx, const Real offset) const noexcept {

test/profile_stellar_collapse.cpp

@@ -63,8 +63,8 @@ struct Bounds {
       }
       min += offset;
       max += offset;
-      min = BDMath::log10(std::abs(min)); // fast, floating-point log
-      max = BDMath::log10(std::abs(max));
+      min = std::log10(std::abs(min)); // fast, floating-point log
+      max = std::log10(std::abs(max));
     }
     grid = RegularGrid1D(min, max, N);
   }

utils/fast-math/logs.hpp

@@ -19,30 +19,52 @@
 #include <cmath>
 #include <ports-of-call/portability.hpp>
 
-// herumi-fmath does not work on device
-// On CPUS it provides another 10% or so speedup
-// for negligible cost in accuracy
-#ifdef PORTABILITY_STRATEGY_KOKKOS
-#define BD_USE_FMATH 0
-#else
-#include <herumi-fmath/fmath.hpp>
-#define BD_USE_FMATH 1
-#endif
 
-namespace BDMath {
 
-  PORTABLE_INLINE_FUNCTION
-  float log10(const float x) {
+namespace singularity {
+namespace Math {
+
+  PORTABLE_FORCEINLINE_FUNCTION
+  double log10(const double x) {
+
     // const double ILOG10 = 1./std::log(10.0);
     constexpr double ILOG10 = 0.43429448190325176;
-    #if BD_USE_FMATH
-    return fmath::log(x)*ILOG10;
-    #else
-    return logf(x)*ILOG10;
-    #endif
+
+    int n{};
+    constexpr const double LOG2OLOG10 = 0.301029995663981195;
+#ifndef SINGULARITY_USE_SINGLE_LOGS
+    const double y = frexp(x, &n); // default is double precision
+#else
+    const float y = frexpf((float)x, &n); // faster but less accurate
+#endif
+#ifdef SINGULARITY_FMATH_USE_ORDER_4
+    // 4th order approximation to log(x) x in [0.5, 1.0)
+    const double expr = -2.36697629372863 + y*(
+			5.2672858163176  + y*(
+                       -5.1242620871906  +   (
+                        2.91607506431871 - 0.69212249971711*y)*y));
+#elif defined(SINGULARITY_FMATH_USE_ORDER_5)
+    // 5th order approximation to log(x) x in [0.5, 1.0)
+    const double expr = -2.5396743497743 + y*(
+                         6.420250933346  + y*(
+                        -8.150589926467  + y*(
+                         6.830564334687  + (
+		        -3.192132453902  + 0.6315814621098*y)*y)));
+#else
+    // 7the order approximation to log(x) x in [0.5, 1.0)
+    const double expr = -2.808428971 + y*(
+                         8.648808309 + y*(
+                       -15.910426569 + y*(
+                        21.53943522  + y*(
+                       -19.56870772  + y*(
+                        11.318698784 +   (
+		        -3.770730277 + 0.5513512194*y)*y)))));
+#endif // SINGULARITY_FMATH_USE_ORDER
+    // log10 x using frexp
+    return ILOG10*expr + n*LOG2OLOG10;
   }
 
-} // BDMath
+} // Math
+} // singularity
 
-#undef BD_USE_FMATH
 #endif //  _SINGULARITY_EOS_UTILS_FAST_MATH_LOGS_