WIP: CREST EoS

singularity-eos/eos/eos.hpp

@@ -33,6 +33,7 @@
 #include <singularity-eos/base/variadic_utils.hpp>
 
 // EOS models
+#include <singularity-eos/eos/eos_crest.hpp>
 #include <singularity-eos/eos/eos_davis.hpp>
 #include <singularity-eos/eos/eos_eospac.hpp>
 #include <singularity-eos/eos/eos_gruneisen.hpp>
@@ -66,7 +67,8 @@ 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
+       SAP_Polynomial,
+       CRESTReactants
 #ifdef SPINER_USE_HDF
        ,
        SpinerEOSDependsRhoT, SpinerEOSDependsRhoSie, StellarCollapse

singularity-eos/eos/eos_crest.hpp

@@ -0,0 +1,204 @@
+//------------------------------------------------------------------------------
+// © 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_CREST_HPP_
+#define _SINGULARITY_EOS_EOS_EOS_CREST_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 CRESTReactants : public EosBase<CRESTReactants> {
+ public:
+  CRESTReactants() = default;
+  PORTABLE_INLINE_FUNCTION
+  CRESTReactants(const Real rho0, const Real a0, const Real a1, const Real a2c,
+                 const Real a2e, const Real a3, const Real b0, const Real b1,
+                 const Real b2c, const Real b2e, const Real b3)
+      : _rho0(rho0), _a0(a0), _a1(a1), _a2c(a2c), _a2e(a2e), _a3(a3), _b0(b0), _b1(b1),
+        _b2c(b2c), _b2e(b2e), _b3(b3) {
+    checkParams();
+  }
+
+  CRESTReactants GetOnDevice() { return *this; }
+  PORTABLE_INLINE_FUNCTION void checkParams() const {
+    PORTABLE_ALWAYS_REQUIRE(_rho0 >= 0, "Reference density must be non-negative");
+  }
+  PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
+    return 0.0;
+  }
+  PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
+    return 0.0;
+  }
+  PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
+    return 0.0;
+  }
+  PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
+    return 0.0;
+  }
+  PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
+    return 0.0;
+  }
+  PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
+    return 0.0;
+  }
+  PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
+    return 0.0;
+  }
+  PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
+    return 0.0;
+  }
+  PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
+    return 0.0;
+  }
+  PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
+    return 0.0;
+  }
+  PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
+    return 0.0;
+  }
+  PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
+    return 0.0;
+  }
+
+  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 = 0.0;
+    sie = 0.0;
+    press = 0.0;
+    cv = 0.0;
+    bmod = 0.0;
+    dpde = 0.0;
+    dvdt = 0.0;
+  }
+  // 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(CRESTReactants)
+  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("CRESTReactants EOS Parameters:\n");
+    printf("    rho0  = %g\n", _rho0);
+    printf("      a0  = %g\n", _a0);
+    printf("      a1  = %g\n", _a1);
+    printf("      a2c = %g\n", _a2c);
+    printf("      a2e = %g\n", _a2e);
+    printf("      a3  = %g\n", _a3);
+    printf("      b0  = %g\n", _b0);
+    printf("      b1  = %g\n", _b1);
+    printf("      b2c = %g\n", _b2c);
+    printf("      b2e = %g\n", _b2e);
+    printf("      b3  = %g\n", _b3);
+  }
+  PORTABLE_INLINE_FUNCTION void
+  DensityEnergyFromPressureTemperature(const Real press, const Real temp, Real *lambda,
+                                       Real &rho, Real &sie) const {
+    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
+    sie = 0.0;
+    rho = 0.0;
+  }
+  inline void Finalize() {}
+  static std::string EosType() { return std::string("CRESTReactants"); }
+  static std::string EosPyType() { return EosType(); }
+
+ private:
+  Real _a0, _a1, _a2c, _a2e, _a3, _b0, _b1, _b2c, _b2e, _b3;
+  // reference values
+  Real _rho0;
+
+  static constexpr const unsigned long _preferred_input =
+      thermalqs::density | thermalqs::specific_internal_energy;
+};
+
+PORTABLE_INLINE_FUNCTION
+void CRESTReactants::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 = 0.0;
+  }
+  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_CREST_HPP_

singularity-eos/eos/singularity_eos.f90

@@ -38,6 +38,7 @@ module singularity_eos
     init_sg_DavisProducts_f,&
     init_sg_DavisReactants_f,&
     init_sg_SAP_Polynomial_f,&
+    init_sg_CRESTReactants_f,&
     init_sg_StiffGas_f,&
     init_sg_SpinerDependsRhoT_f,&
     init_sg_SpinerDependsRhoSie_f,&
@@ -150,6 +151,20 @@ end function init_sg_StiffGas
     end function init_sg_SAP_Polynomial
   end interface
 
+  interface
+    integer(kind=c_int) function &
+      init_sg_CRESTReactants(matindex, eos, rho0, a0, a1, a2c, a2e, a3, b0, b1,&
+                             b2c, b2e, b3, sg_mods_enabled, sg_mods_values) &
+      bind(C, name='init_sg_CRESTReactants')
+      import
+      integer(c_int), value, intent(in)      :: matindex
+      type(c_ptr), value, intent(in)         :: eos
+      real(kind=c_double), value, intent(in) :: rho0, a0, a1, a2c, a2e, a3,&
+                                                b0, b1, b2c, b2e, b3
+      type(c_ptr), value, intent(in)         :: sg_mods_enabled, sg_mods_values
+    end function init_sg_CRESTReactants
+  end interface
+
   interface
     integer(kind=c_int) function &
       init_sg_SpinerDependsRhoT(matindex, eos, filename, id, sg_mods_enabled, &
@@ -388,6 +403,21 @@ integer function init_sg_SAP_Polynomial_f(matindex, eos, rho0, a0, a1, a2c, &
                                  c_loc(sg_mods_enabled), c_loc(sg_mods_values))
   end function init_sg_SAP_Polynomial_f
 
+  integer function init_sg_CRESTReactants_f(matindex, eos, rho0, a0, a1, a2c, &
+                                            a2e, a3, b0, b1, b2c, b2e, b3, &
+                                            sg_mods_enabled, sg_mods_values) &
+    result(err)
+    integer(c_int), value, intent(in) :: matindex
+    type(sg_eos_ary_t), intent(in)    :: eos
+    real(kind=8), value, intent(in)   :: rho0, a0, a1, a2c, a2e, a3,&
+             b0, b1, b2c, b2e, b3
+    integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled
+    real(kind=8), dimension(:), target, intent(inout)        :: sg_mods_values
+    err = init_sg_CRESTReactants(matindex-1, eos%ptr, rho0, a0, a1, a2c, a2e, &
+                                 a3, b0, b1, b2c, b2e, b3, &
+                                 c_loc(sg_mods_enabled), c_loc(sg_mods_values))
+  end function init_sg_CRESTReactants_f
+
   integer function init_sg_StiffGas_f(matindex, eos, gm1, Cv, &
                                       Pinf, qq, &
                                       sg_mods_enabled, sg_mods_values) &

test/CMakeLists.txt

@@ -18,6 +18,7 @@ add_executable(eos_unit_tests
   test_eos_unit.cpp
   test_eos_gruneisen.cpp
   test_eos_sap_polynomial.cpp
+  test_eos_crest.cpp
   test_eos_vinet.cpp
   test_eos_stiff.cpp
 )