lanl · Yurlungur · Aug 25, 2023 · Jun 20, 2023 · Jun 20, 2023 · Jun 20, 2023
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
@@ -35,6 +35,8 @@ jobs:
                   -DSINGULARITY_BUILD_EXAMPLES=ON \
                   -DSINGULARITY_BUILD_PYTHON=ON \
                   -DSINGULARITY_TEST_SESAME=OFF \
+                  -DSINGULARITY_USE_HELMHOLTZ=ON \
+                  -DSINGULARITY_TEST_HELMHOLTZ=ON \
                   -DSINGULARITY_FORCE_SUBMODULE_MODE=ON \
                   ..
                   #-DSINGULARITY_TEST_PYTHON=ON \

diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -23,6 +23,7 @@
 ### Added (new features/APIs/variables/...)
 - [[PR279]](https://github.com/lanl/singularity-eos/pull/279) added noble-abel EoS
 - [[PR274]](https://github.com/lanl/singularity-eos/pull/274) added a stiffened gas EoS
+- [[PR264]](https://github.com/lanl/singularity-eos/pull/274) added a Helmholtz EoS
 - [[PR254]](https://github.com/lanl/singularity-eos/pull/254) added ability to peel off modifiers as needed
 - [[PR250]](https://github.com/lanl/singularity-eos/pull/250) added mass fraction output to stellar collapse eos
 - [[PR226]](https://github.com/lanl/singularity-eos/pull/226) added entropy interpolation to stellar collapse eos

diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -69,6 +69,9 @@ cmake_dependent_option(
 cmake_dependent_option(
   SINGUALRITY_BUILD_STELLARCOLLAPSE2SPINER "Compile stellarcollapse2spiner" ON
   "SINGULARITY_USE_SPINER;SINGULARITY_USE_SPINER_WITH_HDF5" OFF)
+cmake_dependent_option(
+  SINGULARITY_USE_HELMHOLTZ "Include Helmholtz equation of state" OFF
+  "SINGULARITY_USE_SPINER;SINGULARITY_USE_SPINER_WITH_HDF5" OFF)
 
 # testing options
 option(SINGULARITY_BUILD_TESTS "Compile tests" OFF)
@@ -81,6 +84,8 @@ cmake_dependent_option(
   "SINGULARITY_BUILD_TESTS;SINGUALRITY_BUILD_STELLARCOLLAPSE2SPINER" OFF)
 cmake_dependent_option(SINGULARITY_TEST_PYTHON "Test the Python bindings" ON
                        "SINGULARITY_BUILD_TESTS;SINGULARITY_BUILD_PYTHON" OFF)
+cmake_dependent_option(SINGULARITY_TEST_HELMHOLTZ "Test the Helmholtz equation of state" ON
+                       "SINGULARITY_BUILD_TESTS;SINGULARITY_USE_HELMHOLTZ" OFF)
 
 # modify flags options
 option(SINGULARITY_BETTER_DEBUG_FLAGS "Better debug flags for singularity" ON)
@@ -233,6 +238,9 @@ endif()
 if(SINGULARITY_BUILD_CLOSURE)
   target_compile_definitions(singularity-eos PRIVATE SINGULARITY_BUILD_CLOSURE)
 endif()
+if(SINGULARITY_USE_HELMHOLTZ)
+  target_compile_definitions(singularity-eos PUBLIC SINGULARITY_USE_HELMHOLTZ)
+endif()
 
 # ------------------------------------------------------------------------------#
 # Handle dependencies

diff --git a/cmake/install.cmake b/cmake/install.cmake
@@ -91,3 +91,14 @@ export(
   EXPORT singularity-eosTargets
   FILE ${CMAKE_CURRENT_BINARY_DIR}/singularity-eosTargets.cmake
   NAMESPACE singularity-eos::)
+
+# ----------------------------------------------------------------------------#
+# Data files
+# ----------------------------------------------------------------------------#
+
+# install data files needed for various eos models
+if(SINGULARITY_USE_HELMHOLTZ)
+  # install data files
+  install(DIRECTORY ${PROJECT_SOURCE_DIR}/data/helmholtz
+          DESTINATION ${CMAKE_INSTALL_DATADIR}/singularity-eos)
+endif()
diff --git a/data/helmholtz/helm_table.dat b/data/helmholtz/helm_table.dat
diff --git a/doc/sphinx/src/building.rst b/doc/sphinx/src/building.rst
@@ -153,8 +153,22 @@ preconditions:
  ``SINGULARITY_TEST_SESAME``                    ``SINGULARITY_BUILD_TESTS=ON`` ``SINGULARITY_BUILD_SESAME2SPINER=ON``             Test the Sesame table readers.
  ``SINGULARITY_TEST_STELLAR_COLLAPSE``          ``SINGULARITY_BUILD_TESTS=ON`` ``SINGULARITY_BUILD_STELLARCOLLAPSE2SPINER=ON``     Test the Stellar Collapse table readers.
  ``SINGULARITY_TEST_PYTHON``                    ``SINGULARITY_BUILD_TESTS=ON`` ``SINGULARITY_BUILD_PYTHON=ON``                    Test the Python bindings.
+ ``SINGULARITY_USE_HELMHOLTZ``                  ``SINGULARITY_USE_SPINER=ON`` ``SINGULARITY_USE_SPINER_WITH_HDF5=ON``             Use Helmholtz equation of state.
+ ``SINGULARITY_TEST_HELMHOLTZ``                 ``SINGULARITY_USE_HELMHOLTZ``                                                     Build Helmholtz equation of state tests.
 ============================================== ================================================================================= ===========================================
 
+When installing ``singularity-eos``, data files are also installed. The
+follwing options control where the data files are installed:
+
+====================================== ======= ===========================================
+  Option                               Default  Comment
+====================================== ======= ===========================================
+``CMAKE_INSTALL_DATADIR``              <none>  Install directory for data files.
+``CMAKE_INSTALL_DATAROOTDIR``          share   Fallback data install directory.
+====================================== ======= ===========================================
+
+The paths specified by these options are relative to the install prefix.
+
 CMake presets
 -------------
 

diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst
@@ -1288,6 +1288,101 @@ indicates log base 10.
 
 .. _median filter: https://en.wikipedia.org/wiki/Median_filter
 
+
+
+Helmholtz EOS
+``````````````
+
+This is a performance portable implementation of the Helmholtz
+equation of state provided by `Timmes and Swesty`_.  The Helmholtz EOS
+is a three part thermodynamically consistent EOS for a hot, ionized
+gas. It consists of a thermal radiation term:
+
+.. math::
+
+   P = \sigma \cdot T^4
+
+an ions term, treated as an ideal gas:
+
+.. math::
+
+   P = (\gamma - 1) \cdot \rho \cdot e
+
+and a degenerate electron term. Additionally, coulomb force
+corrections can be applied on top of the full model.  This
+multi-component model depends on the relative abundances of electrons
+and ions, as well as the atomic mass and charge of the ions. As such,
+the Helmholtz EOS requires two additional indepenent variables, the
+average atomic mass, Abar, and the average atomic number, Zbar. These
+are passed in through the lambda pointer. As with the other tabulated
+EOS's, the log of the temperature is also stored in the lambda pointer
+as a cache for root finding. Helmholtz provides an enum for indexing
+into the lambda:
+
+* ``Helmholtz::Lambda::Abar`` indexes into the ``Abar`` component of
+  the lambda array.
+* ``Helmholtz::Lambda::Zbar`` indexes into the ``Zbar`` component of
+  the lambda array.
+* ``Helmholtz::Lambda::lT`` indexes into the log temperature cache
+  inside the lambda array.
+
+The degenerate electron term is computed via thermodynamic derivatives
+of the Helmholtz free energy (hence the name Helmholtz EOS). The free
+energy is pre-computed via integrals over the Fermi sphere and
+tabulated in a file provided from `Frank Timmes's website`_.
+
+The table is a simple small ascii file. To ensure thermodyanic
+consistency, the table is interpolated using either biquintic or
+bicubic Hermite polynomials, which are sufficiently high order that
+their high-order derivatives match the underlying data.
+
+.. warning::
+
+   Only a modified version of the table is supported due to the fixed
+   number of colums in the table. This may change in the future.
+   The original table found on `Frank Timmes's website`_ is not supported.
+   A compatible version of the table can be found in the
+   ``data/helmholtz`` directory of the source code, or the data
+   directory specified in the installation configuration.
+
+.. note::
+
+   The implication of interpolating from the free energy is that each
+   EOS evaluation provides ALL relevant EOS data and thermodynamic
+   derivatives. Thus the per-quantity EOS calls are relatively
+   inefficient, and it is instead better to use the FillEos call to
+   get the entire model at once.
+
+The Helmholtz EOS is instantiated by passing in the path to the
+relevant table:
+
+.. code-block::
+
+   Helmholtz(const std::string &filename)
+
+Note that e.g. the Gruneisen parameter is defined differently
+compared to other EOSs. Here the Gruneisen parameter is the
+:math:`\Gamma_3` of Cox & Giuli 1968 - Princiiples of Stellar Structure
+(c&g in the following). Specifically:
+
+.. math::
+
+    \Gamma_3 - 1 = \left. \frac{\mathrm{d} \ln T}{ \mathrm{d} \ln \rho}\right|_\mathrm{ad}
+
+Some important formulas to be used when using this EOS:
+ - the temperature and density exponents (c&g 9.81 9.82)  
+ - the specific heat at constant volume (c&g 9.92)  
+ - the third adiabatic exponent (c&g 9.93)  
+ - the first adiabatic exponent (c&g 9.97)  
+ - the second adiabatic exponent (c&g 9.105)  
+ - the specific heat at constant pressure (c&g 9.98)  
+ - and relativistic formula for the sound speed (c&g 14.29)  
+
+
+.. _Timmes and Swesty: https://doi.org/10.1086/313304
+
+.. _Frank Timmes's website: https://cococubed.com/code_pages/eos.shtml
+
 EOSPAC EOS
 ````````````
 

diff --git a/python/module.cpp b/python/module.cpp
@@ -110,6 +110,11 @@ PYBIND11_MODULE(singularity_eos, m) {
     .def_property_readonly("YeMax", &StellarCollapse::YeMax)
     .def_property_readonly("sieMin", &StellarCollapse::sieMin)
     .def_property_readonly("sieMax", &StellarCollapse::sieMax);
+
+  eos_class<Helmholtz>(m, "Helmholtz")
+    .def(py::init())
+    .def(py::init<std::string&>())
+    .def(py::init<std::string&,bool,bool,bool,bool,bool>());
 #endif
 
 #ifdef SINGULARITY_USE_EOSPAC

diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt
@@ -32,6 +32,7 @@ set(EOS_HEADERS
     eos/eos_vinet.hpp
     eos/eos_builder.hpp
     eos/eos_jwl.hpp
+    eos/eos_helmholtz.hpp
     eos/modifiers/relativistic_eos.hpp
     eos/modifiers/scaled_eos.hpp
     eos/modifiers/ramps_eos.hpp

diff --git a/singularity-eos/base/hermite.hpp b/singularity-eos/base/hermite.hpp
@@ -0,0 +1,99 @@
+//------------------------------------------------------------------------------
+// © 2023. 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.
+//------------------------------------------------------------------------------
+
+#ifndef SINGULARITY_EOS_BASE_HERMITE_HPP_
+#define SINGULARITY_EOS_BASE_HERMITE_HPP_
+
+#include <ports-of-call/portability.hpp>
+#include <singularity-eos/base/math_utils.hpp>
+
+namespace singularity {
+namespace hermite {
+
+// psi0 and its derivatives
+PORTABLE_INLINE_FUNCTION Real psi0(Real z) {
+  return z * z * z * (z * (-6.0 * z + 15.0) - 10.0) + 1.0;
+}
+PORTABLE_INLINE_FUNCTION Real dpsi0(Real z) {
+  return z * z * (z * (-30.0 * z + 60.0) - 30.0);
+}
+PORTABLE_INLINE_FUNCTION Real ddpsi0(Real z) {
+  return z * (z * (-120.0 * z + 180.0) - 60.0);
+}
+
+// psi1 and its derivatives
+PORTABLE_INLINE_FUNCTION Real psi1(Real z) {
+  return z * (z * z * (z * (-3.0 * z + 8.0) - 6.0) + 1.0);
+}
+PORTABLE_INLINE_FUNCTION Real dpsi1(Real z) {
+  return z * z * (z * (-15.0 * z + 32.0) - 18.0) + 1.0;
+}
+PORTABLE_INLINE_FUNCTION Real ddpsi1(Real z) {
+  return z * (z * (-60.0 * z + 96.0) - 36.0);
+}
+
+// psi2 and its derivatives
+PORTABLE_INLINE_FUNCTION Real psi2(Real z) {
+  return 0.5 * z * z * (z * (z * (-z + 3.0) - 3.0) + 1.0);
+}
+PORTABLE_INLINE_FUNCTION Real dpsi2(Real z) {
+  return 0.5 * z * (z * (z * (-5.0 * z + 12.0) - 9.0) + 2.0);
+}
+PORTABLE_INLINE_FUNCTION Real ddpsi2(Real z) {
+  return 0.5 * (z * (z * (-20.0 * z + 36.0) - 18.0) + 2.0);
+}
+
+// biquintic hermite polynomial
+PORTABLE_INLINE_FUNCTION Real h5(Real fi[36], Real w0t, Real w1t, Real w2t, Real w0mt,
+                                 Real w1mt, Real w2mt, Real w0d, Real w1d, Real w2d,
+                                 Real w0md, Real w1md, Real w2md) {
+  return fi[0] * w0d * w0t + fi[1] * w0md * w0t + fi[2] * w0d * w0mt +
+         fi[3] * w0md * w0mt + fi[4] * w0d * w1t + fi[5] * w0md * w1t +
+         fi[6] * w0d * w1mt + fi[7] * w0md * w1mt + fi[8] * w0d * w2t +
+         fi[9] * w0md * w2t + fi[10] * w0d * w2mt + fi[11] * w0md * w2mt +
+         fi[12] * w1d * w0t + fi[13] * w1md * w0t + fi[14] * w1d * w0mt +
+         fi[15] * w1md * w0mt + fi[16] * w2d * w0t + fi[17] * w2md * w0t +
+         fi[18] * w2d * w0mt + fi[19] * w2md * w0mt + fi[20] * w1d * w1t +
+         fi[21] * w1md * w1t + fi[22] * w1d * w1mt + fi[23] * w1md * w1mt +
+         fi[24] * w2d * w1t + fi[25] * w2md * w1t + fi[26] * w2d * w1mt +
+         fi[27] * w2md * w1mt + fi[28] * w1d * w2t + fi[29] * w1md * w2t +
+         fi[30] * w1d * w2mt + fi[31] * w1md * w2mt + fi[32] * w2d * w2t +
+         fi[33] * w2md * w2t + fi[34] * w2d * w2mt + fi[35] * w2md * w2mt;
+}
+
+// cubic hermite polynomial
+// psi0 and its derivatives
+PORTABLE_INLINE_FUNCTION Real xpsi0(Real z) { return z * z * (2.0 * z - 3.0) + 1.0; }
+PORTABLE_INLINE_FUNCTION Real xdpsi0(Real z) { return z * (6.0 * z - 6.0); }
+
+// psi1 and its derivatives
+PORTABLE_INLINE_FUNCTION Real xpsi1(Real z) { return z * (z * (z - 2.0) + 1.0); }
+
+PORTABLE_INLINE_FUNCTION Real xdpsi1(Real z) { return z * (3.0 * z - 4.0) + 1.0; }
+
+// bicubic hermite polynomial
+PORTABLE_INLINE_FUNCTION Real h3(Real fi[16], Real w0t, Real w1t, Real w0mt, Real w1mt,
+                                 Real w0d, Real w1d, Real w0md, Real w1md) {
+  return fi[0] * w0d * w0t + fi[1] * w0md * w0t + fi[2] * w0d * w0mt +
+         fi[3] * w0md * w0mt + fi[4] * w0d * w1t + fi[5] * w0md * w1t +
+         fi[6] * w0d * w1mt + fi[7] * w0md * w1mt + fi[8] * w1d * w0t +
+         fi[9] * w1md * w0t + fi[10] * w1d * w0mt + fi[11] * w1md * w0mt +
+         fi[12] * w1d * w1t + fi[13] * w1md * w1t + fi[14] * w1d * w1mt +
+         fi[15] * w1md * w1mt;
+}
+
+} // namespace hermite
+} // namespace singularity
+
+#endif // SINGULARITY_EOS_BASE_MATH_UTILS_HPP_
diff --git a/singularity-eos/base/math_utils.hpp b/singularity-eos/base/math_utils.hpp
@@ -22,32 +22,40 @@ namespace math_utils {
 
 /*
  * Replaces "SQUARE, CUBE, etc
- * Generic is recursive, but we specialize a few of them strategically
- * A generalization would be to treat even powers specially
+ * Use ternary operator instead of if constexpr to enable compile-time tail-recursion
+ * for various edge cases. Following accelerants are used:
+ * x^{2*i} = x^i x^i
+ * x^i = i x^{i - 1}
+ * x^{-i} = 1/x^i
  */
-template <size_t P>
-PORTABLE_FORCEINLINE_FUNCTION constexpr auto pow(const Real x) {
-  return x * pow<P - 1>(x);
+template <int P>
+PORTABLE_FORCEINLINE_FUNCTION constexpr Real pow(const Real x) {
+  // defining Ppos is important. Otherwise compiler is not smart
+  // enough to terminate recursion.
+  constexpr int Ppos = (P < 0) ? -P : P;
+  return (P >= 0)
+             ? ((Ppos % 2) ? x * pow<Ppos - 1>(x) : pow<Ppos / 2>(x) * pow<Ppos / 2>(x))
+             : 1 / pow<Ppos>(x);
 }
 
 template <>
-PORTABLE_FORCEINLINE_FUNCTION constexpr auto pow<0>(const Real x) {
+PORTABLE_FORCEINLINE_FUNCTION constexpr Real pow<0>(const Real x) {
   return 1.0;
 }
 
 template <>
-PORTABLE_FORCEINLINE_FUNCTION constexpr auto pow<1>(const Real x) {
+PORTABLE_FORCEINLINE_FUNCTION constexpr Real pow<1>(const Real x) {
   return x;
 }
 
-template <>
-PORTABLE_FORCEINLINE_FUNCTION constexpr auto pow<2>(const Real x) {
-  return x * x;
+template <int P>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto ipow10() {
+  return pow<P>(10);
 }
 
-template <>
-PORTABLE_FORCEINLINE_FUNCTION constexpr auto pow<4>(const Real x) {
-  return pow<2>(x) * pow<2>(x);
+PORTABLE_FORCEINLINE_FUNCTION auto pow10(const Real x) {
+  constexpr Real ln10 = 2.30258509299405e+00;
+  return std::exp(ln10 * x);
 }
 
 } // namespace math_utils