Merge branch 'main' into dholladay00/robust_JWL

.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 \

CHANGELOG.md

@@ -25,6 +25,7 @@
 - [[PR265]](https://github.com/lanl/singularity-eos/pull/265) Add missing UnitSystem modifier combinations to variant and EOSPAC
 - [[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
@@ -48,6 +49,7 @@
 - [[PR245]](https://github.com/lanl/singularity-eos/pull/245) Separating get_sg_eos to other files. Build/compilation improvements, warning fixes/suppression.
 
 ### Removed (removing behavior/API/varaibles/...)
+- [[PR293]](https://github.com/lanl/singularity-eos/pull/293) Removing PTofRE function. This will no longer be callable downstream.
 
 ## Release 1.7.0
 Date: 12/14/2022

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

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/data)
+endif()

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
 -------------
 

doc/sphinx/src/contributing.rst

@@ -491,14 +491,13 @@ The CRTP slass structure and static polymorphism
 ````````````````````````````````````````````````
 
 Each of the EOS models in ``singularity-eos`` inherits from a base class in
-order to centralize default functionality and avoid code duplication. The two
-main examples of this are the vector overloads and the ``PTofRE`` scalar lookup
-function. In the vector overloads, a simple for loop is used to iterate over
+order to centralize default functionality and avoid code duplication. The
+main example of this are the vector overloads.
+In the vector overloads, a simple for loop is used to iterate over
 the set of states provided to the function and then call the scalar version on
-each state. The ``PTofRE`` function is designed to provide a common method for
-getting the needed information for a PTE solve from an EOS. Both of these
-features are general to all types of EOS, but are reliant on specific
-implementations of the EOS lookups. In both cases, these functions provide a
+each state. This feature is
+general to all types of EOS, but reliant on specific
+implementations of the EOS lookups. These functions provide a
 default behaviour that we might also want to override for a given equation of
 state.
 

doc/sphinx/src/models.rst

@@ -1288,6 +1288,117 @@ 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
+
+The constructor for the ``Helmholtz`` EOS class looks like
+
+.. code-block:: cpp
+
+  Helmholtz(const std::string &filename, const bool rad = true,
+            const bool gas = true, const bool coul = true,
+            const bool ion = true, const bool ele = true,
+            const bool verbose = false)
+
+where ``filename`` is the file containing the tabulated model. The
+optional arguments ``rad``, ``gas``, ``coul``, ``ion``, and ``ele``
+specify whether to include the radiation, ideal gas, coulomb correction,
+ionization, and electron contributions, respectively. The default is to
+include all terms. The optional argument ``verbose`` specifies whether to print
+out additional information, e.g. when the root find fails to converge.
+
 EOSPAC EOS
 ````````````
 

doc/sphinx/src/using-eos.rst

@@ -547,12 +547,3 @@ root finder to meaningful bound the root search. Similarly,
 
 returns the density at which pressure is minimized for a given
 temperature. This is again useful for root finds.
-
-Finally the method
-
-.. cpp:function:: void PTofRE(Real &rho, Real &sie, Real *lambda, Real &press, Real &temp, Real &dpdr, Real &dpde, Real &dtdr, Real &dtde) const;
-
-returns pressure and temperature, as well as the thermodynamic
-derivatives of pressure and temperature with respect to density and
-specific internal energy, as a function of density and specific
-internal energy.

python/module.cpp

@@ -110,6 +110,13 @@ PYBIND11_MODULE(singularity_eos, m) {
     .def_property_readonly("YeMax", &StellarCollapse::YeMax)
     .def_property_readonly("sieMin", &StellarCollapse::sieMin)
     .def_property_readonly("sieMax", &StellarCollapse::sieMax);
+
+#ifdef SINGULARITY_USE_HELMHOLTZ
+  eos_class<Helmholtz>(m, "Helmholtz")
+    .def(py::init())
+    .def(py::init<std::string&>())
+    .def(py::init<std::string&,bool,bool,bool,bool,bool>());
+#endif
 #endif
 
 #ifdef SINGULARITY_USE_EOSPAC

python/module.hpp

@@ -344,44 +344,6 @@ py::class_<T> eos_class(py::module_ & m, std::string name) {
 
     .def("MinimumDensity", &T::MinimumDensity)
     .def("MinimumTemperature", &T::MinimumTemperature)
-
-    .def("PTofRE", [](const T & self, Real rho, Real sie, py::array_t<Real> lambda) {
-      EOSState s;
-      s.density = rho;
-      s.specific_internal_energy = sie;
-      self.PTofRE(s.density, s.specific_internal_energy, lambda.mutable_data(), s.pressure, s.temperature, s.dpdr, s.dpde, s.dtdr, s.dtde);
-      return s;
-    }, py::arg("rho"), py::arg("sie"), py::arg("lmbda"))
-    .def("PTofRE", [](const T & self, Real rho, Real sie) {
-      EOSState s;
-      s.density = rho;
-      s.specific_internal_energy = sie;
-      self.PTofRE(s.density, s.specific_internal_energy, nullptr, s.pressure, s.temperature, s.dpdr, s.dpde, s.dtdr, s.dtde);
-      return s;
-    }, py::arg("rho"), py::arg("sie"))
-    .def("PTofRE", [](const T & self, py::array_t<Real> rhos, py::array_t<Real> sies, py::array_t<Real> pressures, py::array_t<Real> temperatures, py::array_t<Real> dpdrs, py::array_t<Real> dpdes, py::array_t<Real> dtdrs, py::array_t<Real> dtdes,
-                      const int num, py::array_t<Real> lambdas) {
-      py::buffer_info lambdas_info = lambdas.request();
-      if (lambdas_info.ndim != 2)
-        throw std::runtime_error("lambdas dimension must be 2!");
-
-      if(lambdas_info.shape[1] > 0) {
-        self.PTofRE(rhos.mutable_data(), sies.mutable_data(),
-                    pressures.mutable_data(), temperatures.mutable_data(), dpdrs.mutable_data(),
-                    dpdes.mutable_data(), dtdrs.mutable_data(), dtdes.mutable_data(), num,
-                    LambdaHelper(lambdas));
-      } else {
-        self.PTofRE(rhos.mutable_data(), sies.mutable_data(),
-                    pressures.mutable_data(), temperatures.mutable_data(), dpdrs.mutable_data(),
-                    dpdes.mutable_data(), dtdrs.mutable_data(), dtdes.mutable_data(), num,
-                    NoLambdaHelper());
-      }
-    }, py::arg("rhos"), py::arg("sies"), py::arg("pressures"), py::arg("temperatures"), py::arg("dpdrs"), py::arg("dpdes"), py::arg("dtdrs"), py::arg("dtdes"), py::arg("num"), py::arg("lmbdas"))
-    .def("PTofRE", [](const T & self, py::array_t<Real> rhos, py::array_t<Real> sies, py::array_t<Real> pressures, py::array_t<Real> temperatures, py::array_t<Real> dpdrs, py::array_t<Real> dpdes, py::array_t<Real> dtdrs, py::array_t<Real> dtdes,
-                      const int num) {
-      self.PTofRE(rhos.mutable_data(), sies.mutable_data(), pressures.mutable_data(), temperatures.mutable_data(), dpdrs.mutable_data(), dpdes.mutable_data(), dtdrs.mutable_data(), dtdes.mutable_data(), num, NoLambdaHelper());
-    }, py::arg("rhos"), py::arg("sies"), py::arg("pressures"), py::arg("temperatures"), py::arg("dpdrs"), py::arg("dpdes"), py::arg("dtdrs"), py::arg("dtdes"), py::arg("num"))
-
     .def_property_readonly("nlambda", &T::nlambda)
     .def_property_readonly_static("PreferredInput", [](py::object) { return T::PreferredInput(); })
     .def("PrintParams", &T::PrintParams)

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