Add general entropy functionality to EOS

.gitignore

@@ -8,3 +8,6 @@ utils/variant
 cmake-gen*
 spack-build-*
 spack-configure-*
+._*
+*.sublime*
+

CHANGELOG.md

@@ -2,30 +2,25 @@
 
 ## Current develop
 
-### Fixed (Repair, bugs, ect)
-- [[PR232]](https://github.com/lanl/singularity-eos/pull/228) Fixed uninitialized cmake path variables
-
-### Fixed (issue #227)
-- [[PR228]](https://github.com/lanl/singularity-eos/pull/228) and [[PR229]](https://github.com/lanl/singularity-eos/pull/229) added untracked header files in cmake
-
-### Added (entropy calculation for stellar collapse eos)
-- [[PR226]](https://github.com/lanl/singularity-eos/pull/226) added entropy interpolation to stellar collapse eos
-
 ### Fixed (Repair bugs, etc)
 - [[PR228]](https://github.com/lanl/singularity-eos/pull/228) added untracked header files in cmake
 - [[PR215]](https://github.com/lanl/singularity-eos/pull/215) and [[PR216]](https://github.com/lanl/singularity-eos/pull/216) fix duplicate definition of EPS and fix CI
+- [[PR232]](https://github.com/lanl/singularity-eos/pull/228) Fixed uninitialized cmake path variables
 
 ### Added (new features/APIs/variables/...)
 - [[PR202]](https://github.com/lanl/singularity-eos/pull/202) added the Vinet analytical EOS wth test cases and documentation.
 - [[PR226]](https://github.com/lanl/singularity-eos/pull/226) added entropy interpolation to stellar collapse eos
 - [[PR209]](https://github.com/lanl/singularity-eos/pull/209) added more documentation around how to contribute to the project and also what a contributor can expect from the core team
 - [[PR214]](https://github.com/lanl/singularity-eos/pull/214) added documentation about adding a new EOS
+- [[PR228]](https://github.com/lanl/singularity-eos/pull/228) and [[PR229]](https://github.com/lanl/singularity-eos/pull/229) added untracked header files in cmake
+- [[PR233]](https://github.com/lanl/singularity-eos/pull/233) Added entropy infrastructure to singularity. Add entropy for the ideal gas and modifiers while throwing an error for EOS where entropy is not implemented yet
 - [[PR177]](https://github.com/lanl/singularity-eos/pull/177) added EOSPAC vector functions
 
 ### Changed (changing behavior/API/variables/...)
 - [[PR223]](https://github.com/lanl/singularity-eos/pull/223) Update ports-of-call and add portable error handling
 - [[PR234]](https://github.com/lanl/singularity-eos/pull/234) update ports-of-call to correct for undefined behavior in error handling
 - [[PR219]](https://github.com/lanl/singularity-eos/pull/219) Removed static analysis from re-git pipeline
+- [[PR233]](https://github.com/lanl/singularity-eos/pull/233) Exposed entropy for the EOS type (now required for future EOS)
 
 ### Infrastructure (changes irrelevant to downstream codes)
 - [[PR190]](https://github.com/lanl/singularity-eos/pull/190) update CI on re-git

doc/sphinx/src/modifiers.rst

@@ -65,7 +65,44 @@ energy. The pressure is unchanged under the operation
 
 for some scale parameter :math:`s`.  The ``ScaledEOS`` applies this
 transformation to any equation of state, not just an ideal gas, where
-the pressure may change.
+the pressure may change for different scaling ratios.
+
+Another way of understanding scaling ratios is that the pressure can be written
+as
+
+.. math::
+
+  P = \left(\frac{\partial F}{\partial V} \right)_T
+
+where :math:`F` is the Helmholtz free energy. For a given scaling such that
+:math:`\rho_\mathrm{eos} = s\rho_\mathrm{in}`, the volume obeys the inverse
+scaling. Since the scaling ratio is constant, it can be substituted into the
+above expression so that
+
+.. math::
+
+  P = \left(\frac{\partial F_\mathrm{eos}}{\partial V_\mathrm{eos}} \right)_T
+    = \left(\frac{\partial F_\mathrm{in}}{\partial V_\mathrm{in}} \right)_T
+    = \left(\frac{\partial F_\mathrm{in}}{s \partial V_\mathrm{eos}} \right)_T
+    = \left(\frac{s\partial F_\mathrm{eos}}{s \partial V_\mathrm{eos}} \right)_T
+
+which implies that the Helmholtz free energy must scale in the same way as
+volume (inverse to density) in order to preserve the same pressure. Applying
+this scaling to the definition of the Helmholtz free energy yields
+
+.. math::
+
+  F_\mathrm{eos} = e_\mathrm{eos} - TS_\mathrm{eos} = \frac{1}{R} F_\mathrm{in}
+    = \frac{1}{R}e_\mathrm{in} - T\left(\frac{1}{R}S_\mathrm{in}\right),
+
+where the implicaiton is that this inverse the scaling ratio should also be
+applied to energy. The inverse scaling ratio must be applied to the entropy
+here in order to ensure that all other thermodynamic potentials
+(energy, entropy, and the Gibbs free energy) scale similarly.
+
+where :math:`e` is the internal energy and :math:`S` is the entropy. The
+implication is that the same scaling should be applied to the energy and entropy
+to maintain thermodynamic consistency.
 
 The constructor for ``ScaledEOS`` accepts the underlying model, and
 the scale parameter. For example, a shifted ideal gas might be

doc/sphinx/src/sphinx-doc.rst

@@ -33,22 +33,22 @@ with shared machines or python distributions.
 First, ensure that you are running a modern version of python (i.e. python 3 of
 some flavor)
 
-.. code-block::
-   :language:bash
+.. code-block:: bash
+
    $ python --version
    Python 3.9.7
 
 Then, use pip to install :code:`spinx` and the RTD theme
 
-.. code-block::
-   :language:bash
+.. code-block:: bash
 
    pip install --user sphinx sphinx-rtd-theme
 
 Now, navigate to the :code:`../doc/sphinx` directory where a :code:`make help`
 shows all of the available ways to build the documentation
 
 .. code-block::
+
    :language:bash
 
    $ make help

singularity-eos/eos/eos_base.hpp

@@ -15,9 +15,11 @@
 #ifndef _SINGULARITY_EOS_EOS_EOS_BASE_
 #define _SINGULARITY_EOS_EOS_EOS_BASE_
 
+#include <cstring>
 #include <string>
 
 #include <ports-of-call/portability.hpp>
+#include <ports-of-call/portable_errors.hpp>
 
 namespace singularity {
 namespace mfuncname {
@@ -50,7 +52,10 @@ namespace eos_base {
   using EosBase<EOSDERIVED>::MinimumDensity;                                             \
   using EosBase<EOSDERIVED>::MinimumTemperature;                                         \
   using EosBase<EOSDERIVED>::PTofRE;                                                     \
-  using EosBase<EOSDERIVED>::FillEos;
+  using EosBase<EOSDERIVED>::FillEos;                                                    \
+  using EosBase<EOSDERIVED>::EntropyFromDensityTemperature;                              \
+  using EosBase<EOSDERIVED>::EntropyFromDensityInternalEnergy;                           \
+  using EosBase<EOSDERIVED>::EntropyIsNotEnabled;
 
 /*
 This is a CRTP that allows for static inheritance so that default behavior for
@@ -150,6 +155,48 @@ class EosBase {
     PressureFromDensityInternalEnergy(rhos, sies, pressures, num, lambdas);
   }
   template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void EntropyFromDensityTemperature(ConstRealIndexer &&rhos,
+                                            ConstRealIndexer &&temperatures,
+                                            RealIndexer &&entropies, const int num,
+                                            LambdaIndexer &&lambdas) const {
+    static auto const name = SG_MEMBER_FUNC_NAME();
+    static auto const cname = name.c_str();
+    CRTP copy = *(static_cast<CRTP const *>(this));
+    portableFor(
+        cname, 0, num, PORTABLE_LAMBDA(const int i) {
+          entropies[i] =
+              copy.EntropyFromDensityTemperature(rhos[i], temperatures[i], lambdas[i]);
+        });
+  }
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void
+  EntropyFromDensityTemperature(ConstRealIndexer &&rhos, ConstRealIndexer &&temperatures,
+                                RealIndexer &&entropies, Real * /*scratch*/,
+                                const int num, LambdaIndexer &&lambdas) const {
+    EntropyFromDensityTemperature(rhos, temperatures, entropies, num, lambdas);
+  }
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void EntropyFromDensityInternalEnergy(ConstRealIndexer &&rhos,
+                                               ConstRealIndexer &&sies,
+                                               RealIndexer &&entropies, const int num,
+                                               LambdaIndexer &&lambdas) const {
+    static auto const name = SG_MEMBER_FUNC_NAME();
+    static auto const cname = name.c_str();
+    CRTP copy = *(static_cast<CRTP const *>(this));
+    portableFor(
+        cname, 0, num, PORTABLE_LAMBDA(const int i) {
+          entropies[i] =
+              copy.EntropyFromDensityInternalEnergy(rhos[i], sies[i], lambdas[i]);
+        });
+  }
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void
+  EntropyFromDensityInternalEnergy(ConstRealIndexer &&rhos, ConstRealIndexer &&sies,
+                                   RealIndexer &&entropies, Real * /*scratch*/,
+                                   const int num, LambdaIndexer &&lambdas) const {
+    EntropyFromDensityInternalEnergy(rhos, sies, entropies, num, lambdas);
+  }
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
   inline void SpecificHeatFromDensityTemperature(ConstRealIndexer &&rhos,
                                                  ConstRealIndexer &&temperatures,
                                                  RealIndexer &&cvs, const int num,
@@ -339,69 +386,15 @@ class EosBase {
   PORTABLE_INLINE_FUNCTION
   Real RhoPmin(const Real temp) const { return 0.0; }
 
-  // Specialzied vector version of PTofRE maybe more suited for EOSPAC
-  // }
-  // template<typename RealIndexer, typename LambdaIndexer>
-  // inline
-  // void PTofRE(RealIndexer &&rhos, RealIndexer &&sies,
-  //             RealIndexer &&presses, RealIndexer &&temps,
-  //             RealIndexer &&dpdrs, RealIndexer &&dpdes,
-  //             RealIndexer &&dtdrs, RealIndexer &&dtdes,
-  //             const int num, LambdaIndexer &&lambdas) const {
-  //   // Get the dervived class
-  //   auto eos = static_cast<CRTP const&>(*this);
-
-  //   // Lookup at density and energy
-  //   eos.PressureFromDensityInternalEnergy(rhos, sies, num, presses,
-  //                                               lambdas);
-  //   eos.TemperatureFromDensityInternalEnergy(rhos, sies, num, temps,
-  //                                                  lambdas);
-  //   // Peturbation factor
-  //   Real factor = 1. + 1.0e-06;
-
-  //   // Perturb densities and do lookups
-  //   portableFor(
-  //       'PerturbDensities', 0, num, PORTABLE_LAMBDA(const int i) {
-  //         rhos[i] *= factor;
-  //       }
-  //   );
-  //   eos.PressureFromDensityInternalEnergy(rhos, sies, num, dpdrs,
-  //                                               lambdas);
-  //   eos.TemperatureFromDensityInternalEnergy(rhos, sies, num, dtdrs,
-  //                                                  lambdas);
-
-  //   // Reset densities, perturb energies, and do lookups
-  //   portableFor(
-  //       'PerturbEnergiesResetDensities', 0, num,
-  //       PORTABLE_LAMBDA(const int i) {
-  //         sies[i] *= factor;
-  //         rhos[i] /= factor;
-  //       }
-  //   );
-  //   eos.PressureFromDensityInternalEnergy(rhos, sies, dpdes, num,
-  //                                               lambdas);
-  //   eos.TemperatureFromDensityInternalEnergy(rhos, sies, dtdes, num,
-  //                                                  lambdas);
-
-  //   // Reset the energies to their original values
-  //   portableFor(
-  //       'ResetEnergies', 0, num, PORTABLE_LAMBDA(const int i) {
-  //         sies[i] /= factor;
-  //       }
-  //   );
-
-  //   // Calculate the derivatives
-  //   portableFor(
-  //       'CalculateDerivatives', 0., num, PORTABLE_LAMBDA(const int i) {
-  //         dpdrs[i] = (dpdrs[i] - presses[i]) / (rhos[i] * (1. - factor));
-  //         dpdes[i] = (dpdes[i] - presses[i]) / (sies[i] * (1. - factor));
-  //         dtdrs[i] = (dtdrs[i] - temps[i]) / (rhos[i] * (1. - factor));
-  //         dtdes[i] = (dtdes[i] - temps[i]) / (sies[i] * (1. - factor));
-  //       }
-  //   );
-
-  //   return;
-  // }
+  // Default entropy behavior is to return an error
+  PORTABLE_FORCEINLINE_FUNCTION
+  void EntropyIsNotEnabled(const char *eosname) const {
+    // Construct the error message using char* so it works on device
+    char msg[120] = "Entropy is not enabled for the '";
+    strcat(msg, eosname);
+    strcat(msg, "' EOS");
+    PORTABLE_ALWAYS_THROW_OR_ABORT(msg);
+  }
 };
 } // namespace eos_base
 } // namespace singularity

singularity-eos/eos/eos_davis.hpp

@@ -60,6 +60,16 @@ class DavisReactants : public EosBase<DavisReactants> {
       const Real rho, const Real sie, Real *lambda = nullptr) const {
     return Ps(rho) + Gamma(rho) * rho * (sie - Es(rho));
   }
+  PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    EntropyIsNotEnabled("DavisReactants");
+    return 1.0;
+  }
+  PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    EntropyIsNotEnabled("DavisReactants");
+    return 1.0;
+  }
   PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(
       const Real rho, const Real temp, Real *lambda = nullptr) const {
     return SpecificHeatFromDensityInternalEnergy(
@@ -152,6 +162,16 @@ class DavisProducts : public EosBase<DavisProducts> {
       const Real rho, const Real sie, Real *lambda = nullptr) const {
     return Ps(rho) + rho * Gamma(rho) * (sie - Es(rho));
   }
+  PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    EntropyIsNotEnabled("DavisProducts");
+    return 1.0;
+  }
+  PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    EntropyIsNotEnabled("DavisProducts");
+    return 1.0;
+  }
   PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(
       const Real rho, const Real temperature, Real *lambda = nullptr) const {
     return _Cv;

singularity-eos/eos/eos_eospac.hpp

@@ -54,6 +54,10 @@ class EOSPAC : public EosBase<EOSPAC> {
       const Real rho, const Real temperature, Real *lambda = nullptr) const;
   PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy(
       const Real rho, const Real sie, Real *lambda = nullptr) const;
+  PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const;
+  PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const;
   PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(
       const Real rho, const Real temperature, Real *lambda = nullptr) const;
   PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy(
@@ -82,6 +86,8 @@ class EOSPAC : public EosBase<EOSPAC> {
   using EosBase<EOSPAC>::InternalEnergyFromDensityTemperature;
   using EosBase<EOSPAC>::PressureFromDensityTemperature;
   using EosBase<EOSPAC>::PressureFromDensityInternalEnergy;
+  using EosBase<EOSPAC>::EntropyFromDensityTemperature;
+  using EosBase<EOSPAC>::EntropyFromDensityInternalEnergy;
   using EosBase<EOSPAC>::SpecificHeatFromDensityTemperature;
   using EosBase<EOSPAC>::SpecificHeatFromDensityInternalEnergy;
   using EosBase<EOSPAC>::BulkModulusFromDensityTemperature;
@@ -90,6 +96,7 @@ class EOSPAC : public EosBase<EOSPAC> {
   using EosBase<EOSPAC>::GruneisenParamFromDensityInternalEnergy;
   using EosBase<EOSPAC>::PTofRE;
   using EosBase<EOSPAC>::FillEos;
+  using EosBase<EOSPAC>::EntropyIsNotEnabled;
 
   // EOSPAC vector implementations
   template <typename LambdaIndexer>
@@ -810,6 +817,12 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::PressureFromDensityTemperature(const Real
   return Real(pressureFromSesame(P[0]));
 }
 
+PORTABLE_INLINE_FUNCTION Real EOSPAC::EntropyFromDensityTemperature(
+    const Real rho, const Real temperature, Real *lambda) const {
+  EntropyIsNotEnabled("EOSPAC");
+  return 1.0;
+}
+
 PORTABLE_INLINE_FUNCTION void EOSPAC::FillEos(Real &rho, Real &temp, Real &sie,
                                               Real &press, Real &cv, Real &bmod,
                                               const unsigned long output,
@@ -923,6 +936,12 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::PressureFromDensityInternalEnergy(
   Real temp = TemperatureFromDensityInternalEnergy(rho, sie, lambda);
   return PressureFromDensityTemperature(rho, temp, lambda);
 }
+PORTABLE_INLINE_FUNCTION Real EOSPAC::EntropyFromDensityInternalEnergy(
+    const Real rho, const Real sie, Real *lambda) const {
+  using namespace EospacWrapper;
+  const Real temp = TemperatureFromDensityInternalEnergy(rho, sie, lambda);
+  return EntropyFromDensityTemperature(rho, temp, lambda);
+}
 PORTABLE_INLINE_FUNCTION Real EOSPAC::SpecificHeatFromDensityInternalEnergy(
     const Real rho, const Real sie, Real *lambda) const {
   using namespace EospacWrapper;

singularity-eos/eos/eos_gruneisen.hpp

@@ -76,6 +76,10 @@ class Gruneisen : public EosBase<Gruneisen> {
       const Real rho, const Real temperature, Real *lambda = nullptr) const;
   PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy(
       const Real rho, const Real sie, Real *lambda = nullptr) const;
+  PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const;
+  PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const;
   PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(
       const Real rho, const Real temperatummmmmmre, Real *lambda = nullptr) const {
     return _Cv;
@@ -310,6 +314,12 @@ PORTABLE_INLINE_FUNCTION Real Gruneisen::PressureFromDensityInternalEnergy(
   }
   return P_H + Gamma(rho) * rho * (sie - E_H);
 }
+PORTABLE_INLINE_FUNCTION Real Gruneisen::EntropyFromDensityInternalEnergy(
+    const Real rho_in, const Real sie, Real *lambda) const {
+  const Real rho = std::min(rho_in, _rho_max);
+  EntropyIsNotEnabled("Gruneisen");
+  return 1.0;
+}
 PORTABLE_INLINE_FUNCTION Real Gruneisen::BulkModulusFromDensityInternalEnergy(
     const Real rho_in, const Real sie, Real *lambda) const {
   using namespace gruneisen_utils;
@@ -332,6 +342,12 @@ PORTABLE_INLINE_FUNCTION Real Gruneisen::PressureFromDensityTemperature(
   return PressureFromDensityInternalEnergy(
       rho, InternalEnergyFromDensityTemperature(rho, temp));
 }
+PORTABLE_INLINE_FUNCTION Real Gruneisen::EntropyFromDensityTemperature(
+    const Real rho_in, const Real temp, Real *lambda) const {
+  const Real rho = std::min(rho_in, _rho_max);
+  const Real sie = InternalEnergyFromDensityTemperature(rho, temp);
+  return EntropyFromDensityInternalEnergy(rho, sie);
+}
 PORTABLE_INLINE_FUNCTION Real Gruneisen::BulkModulusFromDensityTemperature(
     const Real rho_in, const Real temp, Real *lambda) const {
   const Real rho = std::min(rho_in, _rho_max);