Bjm/carnstar

CHANGELOG.md

@@ -4,6 +4,7 @@
 - [[PR269]](https://github.com/lanl/singularity-eos/pull/269) Add SAP Polynomial EoS
 
 ### Fixed (Repair bugs, etc)
+- [[PR292]](https://github.com/lanl/singularity-eos/pull/292) Added Carnahan-Starling EoS
 - [[PR290]](https://github.com/lanl/singularity-eos/pull/290) Added target guards on export config
 - [[PR288]](https://github.com/lanl/singularity-eos/pull/288) Don't build tests that depend on spiner when spiner is disabled
 - [[PR287]](https://github.com/lanl/singularity-eos/pull/287) Fix testing logic with new HDF5 options

doc/sphinx/src/models.rst

@@ -542,6 +542,77 @@ these values are not set, they will be the same as those returned by the
 conditions are given, the return values of the :code:`ValuesAtReferenceState()`
 function will not be the same.
 
+Carnahan-Starling
+`````````````````
+
+The (quasi-exact) Carnahan-Starling model in ``singularity-eos`` takes
+the form
+
+.. math::
+
+    P = Z(\rho) \rho (e-q) (\gamma-1)
+
+.. math::
+
+    Z(\rho) = \frac{1+\eta+\eta^2-\eta^3}{(1-\eta)^3}
+
+.. math::
+
+    \eta = b\rho
+
+.. math::
+
+    e = C_V T + q,
+
+the inputs for the equation of state are similar to the Noble-Abel model.
+As with the Noble-Abel EOS, it should be noted that covolume is physically
+significant as it represents the maximum compressibility of the gas,
+and as a result it should be non-negative.
+
+The entropy for the Carnahan-Starling EOS is given by
+
+.. math::
+
+    S =  C_V \ln\left(\frac{T}{T_0}\right) + C_V (\gamma-1) \left\{ \ln\left(\frac{v}
+     {v_0}\right) - S^{CS} \right\} + q',
+
+.. math::
+   S^{CS} = b\left(4\left(\frac{1}{v-b} - \frac{1}{v_0-b}\right)+
+     b\left(\frac{1}{(v-b)^2} - \frac{1}{(v_0-b)^2}\right)\right)
+
+where :math:`S(\rho_0,T_0)=q'`. By default, :math:`T_0 = 298` K and the
+reference density is given by
+
+.. math::
+
+    P_0 = \rho_0 Z(\rho_0) C_V T_0(\gamma-1),
+
+where :math:`P_0` is by default 1 bar. Denisty is obtained via Newton iteration.
+
+The settable parameters for this EOS are :math:`\gamma-1`, specific
+heat capacity (:math:`C_V`), covolume (:math:`b`) and offset internal energy (:math:`q`). Optionally, the reference state for the entropy calculation can
+be provided by setting the reference temperature, pressure, and entropy offset.
+
+The ``CarnahanStarling`` EOS constructor has four arguments: ``gm1``, which is :math:`\gamma-1`; ``Cv``, the
+specific heat :math:`C_V`; :math:`b`, the covolume; and :math:`q`, the internal energy offset.
+
+.. code-block:: cpp
+
+    CarnahanStarling(Real gm1, Real Cv, Real b, Real q)
+
+Optionally, the reference state for the entropy calculation,
+can be provided in the constructor via ``qp``, ``T0`` and ``P0``:
+
+.. code-block:: cpp
+
+    CarnahanStarling(Real gm1, Real Cv, Real b, Real q, Real qp, Real T0, Real P0)
+
+Note that these parameters are provided solely for the entropy calculation. When
+these values are not set, they will be the same as those returned by the
+:code:`ValuesAtReferenceState()` function. However, if the entropy reference
+conditions are given, the return values of the :code:`ValuesAtReferenceState()`
+function will not be the same.
+
 Gruneisen EOS
 `````````````
 

singularity-eos/CMakeLists.txt

@@ -42,6 +42,7 @@ set(EOS_HEADERS
     eos/eos_eospac.hpp
     eos/eos_noble_abel.hpp
     eos/eos_stiff.hpp
+    eos/eos_carnahan_starling.hpp
 )
 
 set(EOS_SRCS

singularity-eos/eos/eos.hpp

@@ -33,6 +33,7 @@
 #include <singularity-eos/base/variadic_utils.hpp>
 
 // EOS models
+#include <singularity-eos/eos/eos_carnahan_starling.hpp>
 #include <singularity-eos/eos/eos_davis.hpp>
 #include <singularity-eos/eos/eos_eospac.hpp>
 #include <singularity-eos/eos/eos_gruneisen.hpp>
@@ -67,7 +68,7 @@ using singularity::detail::transform_variadic_list;
 // all eos's
 static constexpr const auto full_eos_list =
     tl<IdealGas, Gruneisen, Vinet, JWL, DavisReactants, DavisProducts, StiffGas,
-       SAP_Polynomial, NobleAbel
+       SAP_Polynomial, NobleAbel, CarnahanStarling
 #ifdef SINGULARITY_USE_SPINER_WITH_HDF5
        ,
        SpinerEOSDependsRhoT, SpinerEOSDependsRhoSie, StellarCollapse

singularity-eos/eos/eos_carnahan_starling.hpp

@@ -0,0 +1,262 @@
+//------------------------------------------------------------------------------
+// © 2021-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_EOS_EOS_CARNAHAN_STARLING_HPP_
+#define _SINGULARITY_EOS_EOS_EOS_CARNAHAN_STARLING_HPP_
+
+// stdlib
+#include <cmath>
+#include <cstdio>
+#include <string>
+
+// Ports-of-call
+#include <ports-of-call/portability.hpp>
+
+// Base stuff
+#include <singularity-eos/base/constants.hpp>
+#include <singularity-eos/base/eos_error.hpp>
+#include <singularity-eos/base/robust_utils.hpp>
+#include <singularity-eos/eos/eos_base.hpp>
+
+namespace singularity {
+
+using namespace eos_base;
+
+class CarnahanStarling : public EosBase<CarnahanStarling> {
+ public:
+  CarnahanStarling() = default;
+  PORTABLE_INLINE_FUNCTION CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq)
+      : _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(ROOM_TEMPERATURE),
+        _P0(ATMOSPHERIC_PRESSURE), _qp(0.0),
+        _rho0(DensityFromPressureTemperature(_P0, _T0)), _vol0(robust::ratio(1.0, _rho0)),
+        _sie0(Cv * _T0 + qq),
+        _bmod0(_rho0 * Cv * _T0 *
+               (PartialRhoZedFromDensity(_rho0) +
+                ZedFromDensity(_rho0) * ZedFromDensity(_rho0) * gm1)),
+        _dpde0(_rho0 * ZedFromDensity(_rho0) * gm1) {
+    checkParams();
+  }
+  PORTABLE_INLINE_FUNCTION CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq, Real qp,
+                                            Real T0, Real P0)
+      : _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(T0), _P0(P0), _qp(qp),
+        _rho0(DensityFromPressureTemperature(P0, T0)), _vol0(robust::ratio(1.0, _rho0)),
+        _sie0(Cv * T0 + qq),
+        _bmod0(_rho0 * Cv * T0 *
+               (PartialRhoZedFromDensity(_rho0) +
+                ZedFromDensity(_rho0) * ZedFromDensity(_rho0) * gm1)),
+        _dpde0(_rho0 * ZedFromDensity(_rho0) * gm1) {
+    checkParams();
+  }
+  CarnahanStarling GetOnDevice() { return *this; }
+  PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    return std::max(robust::SMALL(), (sie - _qq) / _Cv);
+  }
+  PORTABLE_INLINE_FUNCTION Real ZedFromDensity(const Real rho,
+                                               Real *lambda = nullptr) const {
+    const Real eta = _bb * rho;
+    const Real zed = robust::ratio(1.0 + eta + eta * eta - eta * eta * eta,
+                                   (1.0 - eta) * (1.0 - eta) * (1.0 - eta));
+    return zed;
+  }
+  PORTABLE_INLINE_FUNCTION Real PartialRhoZedFromDensity(const Real rho,
+                                                         Real *lambda = nullptr) const {
+    const Real eta = _bb * rho;
+    return robust::ratio(eta * eta * eta * eta - 4.0 * eta * eta * eta + 4.0 * eta * eta +
+                             4.0 * eta + 1.0,
+                         (1.0 - eta) * (1.0 - eta) * (1.0 - eta) * (1.0 - eta));
+  }
+  PORTABLE_INLINE_FUNCTION Real DensityFromPressureTemperature(
+      const Real press, const Real temperature, Real *lambda = nullptr) const {
+    Real real_root;
+    // iteration related variables
+    int newton_counter = 0;
+    const int newton_max = 50;
+    const Real rel_tol = 1.e-12;
+    real_root = 0.0; // initial guess
+    for (int i = 0; i < newton_max; i++) {
+      Real real_root_old = real_root;
+      Real newton_f =
+          _Cv * temperature * _gm1 * real_root_old * ZedFromDensity(real_root_old) -
+          press;
+      Real newton_f_prime =
+          _Cv * temperature * _gm1 * PartialRhoZedFromDensity(real_root_old);
+      real_root = real_root_old - robust::ratio(newton_f, newton_f_prime);
+      Real rel_err = std::abs(robust::ratio(real_root - real_root_old, real_root));
+      newton_counter++;
+      if (rel_err < rel_tol) {
+        return std::max(robust::SMALL(), real_root);
+      }
+    }
+    if (newton_counter == newton_max) {
+      printf("*** (Warning) DensityFromPressureTemperature :: Convergence not met in "
+             "Carnahan-Starling EoS ***");
+    }
+    return std::max(robust::SMALL(), real_root);
+  }
+  PORTABLE_INLINE_FUNCTION void checkParams() const {
+    PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Heat capacity must be positive");
+    PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive");
+    PORTABLE_ALWAYS_REQUIRE(_bb >= 0, "Covolume must be positive");
+  }
+  PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    return std::max(_qq, _Cv * temperature + _qq);
+  }
+  PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    const Real zed = ZedFromDensity(rho);
+    return std::max(robust::SMALL(), zed * rho * temperature * _Cv * _gm1);
+  }
+  PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    const Real zed = ZedFromDensity(rho);
+    return std::max(robust::SMALL(), zed * rho * (sie - _qq) * _gm1);
+  }
+  PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    const Real vol = robust::ratio(1.0, rho);
+    return _Cv * std::log(robust::ratio(temperature, _T0) + robust::SMALL()) +
+           _gm1 * _Cv * std::log(robust::ratio(vol, _vol0) + robust::SMALL()) -
+           _gm1 * _Cv * _bb *
+               (4.0 * (robust::ratio(1.0, vol - _bb) - robust::ratio(1.0, _vol0 - _bb)) +
+                _bb * (robust::ratio(1.0, (vol - _bb) * (vol - _bb)) -
+                       robust::ratio(1.0, (_vol0 - _bb) * (_vol0 - _bb)))) +
+           _qp;
+  }
+  PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    const Real vol = robust::ratio(1.0, rho);
+    return _Cv * std::log(robust::ratio(sie - _qq, _sie0 - _qq) + robust::SMALL()) +
+           _gm1 * _Cv * std::log(robust::ratio(vol, _vol0) + robust::SMALL()) -
+           _gm1 * _Cv * _bb *
+               (4.0 * (robust::ratio(1.0, vol - _bb) - robust::ratio(1.0, _vol0 - _bb)) +
+                _bb * (robust::ratio(1.0, (vol - _bb) * (vol - _bb)) -
+                       robust::ratio(1.0, (_vol0 - _bb) * (_vol0 - _bb)))) +
+           _qp;
+  }
+  PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    return _Cv;
+  }
+  PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    return _Cv;
+  }
+  PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    return std::max(robust::SMALL(),
+                    rho * _Cv * temperature * _gm1 *
+                        (PartialRhoZedFromDensity(rho) +
+                         ZedFromDensity(rho) * ZedFromDensity(rho) * _gm1));
+  }
+  PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    return std::max(robust::SMALL(),
+                    rho * (sie - _qq) * _gm1 *
+                        (PartialRhoZedFromDensity(rho) +
+                         ZedFromDensity(rho) * ZedFromDensity(rho) * _gm1));
+  }
+  PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    return ZedFromDensity(rho) * _gm1;
+  }
+  PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    return ZedFromDensity(rho) * _gm1;
+  }
+  PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
+                                        Real &cv, Real &bmod, const unsigned long output,
+                                        Real *lambda = nullptr) const;
+  PORTABLE_INLINE_FUNCTION
+  void ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,
+                              Real &bmod, Real &dpde, Real &dvdt,
+                              Real *lambda = nullptr) const {
+    // use STP: 1 atmosphere, room temperature
+    rho = _rho0;
+    temp = _T0;
+    sie = _sie0;
+    press = _P0;
+    cv = _Cv;
+    bmod = _bmod0;
+    dpde = _dpde0;
+  }
+  // Generic functions provided by the base class. These contain e.g. the vector
+  // overloads that use the scalar versions declared here
+  SG_ADD_BASE_CLASS_USINGS(CarnahanStarling)
+  PORTABLE_INLINE_FUNCTION
+  int nlambda() const noexcept { return 0; }
+  static constexpr unsigned long PreferredInput() { return _preferred_input; }
+  static inline unsigned long scratch_size(std::string method, unsigned int nelements) {
+    return 0;
+  }
+  static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; }
+  PORTABLE_INLINE_FUNCTION void PrintParams() const {
+    printf("Carnahan-Starling Parameters:\nGamma = %g\nCv    = %g\nb     = %g\nq     = "
+           "%g\n",
+           _gm1 + 1.0, _Cv, _bb, _qq);
+  }
+  PORTABLE_INLINE_FUNCTION void
+  DensityEnergyFromPressureTemperature(const Real press, const Real temp, Real *lambda,
+                                       Real &rho, Real &sie) const {
+    sie = std::max(_qq, _Cv * temp + _qq);
+    rho = DensityFromPressureTemperature(press, temp);
+  }
+  inline void Finalize() {}
+  static std::string EosType() { return std::string("CarnahanStarling"); }
+  static std::string EosPyType() { return EosType(); }
+
+ private:
+  Real _Cv, _gm1, _bb, _qq;
+  // optional reference state variables
+  Real _T0, _P0, _qp;
+  // reference values
+  Real _rho0, _vol0, _sie0, _bmod0, _dpde0;
+  // static constexpr const Real _T0 = ROOM_TEMPERATURE;
+  // static constexpr const Real _P0 = ATMOSPHERIC_PRESSURE;
+  static constexpr const unsigned long _preferred_input =
+      thermalqs::density | thermalqs::specific_internal_energy;
+};
+
+PORTABLE_INLINE_FUNCTION
+void CarnahanStarling::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,
+                               Real &bmod, const unsigned long output,
+                               Real *lambda) const {
+  if (output & thermalqs::density && output & thermalqs::specific_internal_energy) {
+    if (output & thermalqs::pressure || output & thermalqs::temperature) {
+      UNDEFINED_ERROR;
+    }
+    DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie);
+  }
+  if (output & thermalqs::pressure && output & thermalqs::specific_internal_energy) {
+    if (output & thermalqs::density || output & thermalqs::temperature) {
+      UNDEFINED_ERROR;
+    }
+    sie = InternalEnergyFromDensityTemperature(rho, temp, lambda);
+  }
+  if (output & thermalqs::temperature && output & thermalqs::specific_internal_energy) {
+    sie = robust::ratio(press, ZedFromDensity(rho) * rho * _gm1) + _qq;
+  }
+  if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie);
+  if (output & thermalqs::temperature)
+    temp = TemperatureFromDensityInternalEnergy(rho, sie);
+  if (output & thermalqs::bulk_modulus)
+    bmod = BulkModulusFromDensityInternalEnergy(rho, sie);
+  if (output & thermalqs::specific_heat)
+    cv = SpecificHeatFromDensityInternalEnergy(rho, sie);
+}
+
+} // namespace singularity
+
+#endif // _SINGULARITY_EOS_EOS_EOS_CARNAHAN_STARLING_HPP_

singularity-eos/eos/singularity_eos.cpp

@@ -212,6 +212,25 @@ int init_sg_NobleAbel(const int matindex, EOS *eos, const double gm1, const doub
   return init_sg_NobleAbel(matindex, eos, gm1, Cv, bb, qq, def_en, def_v);
 }
 
+int init_sg_CarnahanStarling(const int matindex, EOS *eos, const double gm1,
+                             const double Cv, const double bb, const double qq,
+                             int const *const enabled, double *const vals) {
+  assert(matindex >= 0);
+  EOS eosi = SGAPPLYMODSIMPLE(CarnahanStarling(gm1, Cv, bb, qq));
+  if (enabled[3] == 1) {
+    singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2],
+                                           vals[3], vals[4], vals[5]);
+  }
+  EOS eos_ = SGAPPLYMOD(CarnahanStarling(gm1, Cv, bb, qq));
+  eos[matindex] = eos_.GetOnDevice();
+  return 0;
+}
+
+int init_sg_CarnahanStarling(const int matindex, EOS *eos, const double gm1,
+                             const double Cv, const double bb, const double qq) {
+  return init_sg_CarnahanStarling(matindex, eos, gm1, Cv, bb, qq, def_en, def_v);
+}
+
 #ifdef SINGULARITY_USE_SPINER_WITH_HDF5
 
 int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename,