Bjm/stiffgas

.github/workflows/deps.yml

@@ -39,7 +39,7 @@ jobs:
           run: git diff --exit-code --compact-summary
           id: no_change
         - name: create branch
-          if: ${{ steps.no_change.conclusion == 'failure' }}
+          if: ${{ failure() && steps.no_change.conclusion == 'failure' }}
           run: |
             git config user.name 'github-actions[bot]'
             git config user.email 'github-actions[bot]@users.noreply.github.com'
@@ -49,7 +49,7 @@ jobs:
             git commit -m "spack updates"
             git push -f --set-upstream origin github-bot/update_spackages
         - name: create pull request
-          if: ${{ steps.no_change.conclusion == 'failure' }}
+          if: ${{ failure() && steps.no_change.conclusion == 'failure' }}
           run: gh pr create -B main -H github-bot/update_spackages --title 'Update spackages' --body 'Created by Github action'
           env:
             GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}

.gitlab-ci.yml

@@ -20,6 +20,7 @@ variables:
   # DEVEL_SPACK_CHECKOUT: "${CI_PROJECT_DIR}/spack"
 
 before_script:
+  - export PYTHONNOUSERSITE=1
   - export SPACK_DISABLE_LOCAL_CONFIG=true
   - export SPACK_USER_CACHE_PATH=${SPACK_USER_CACHE_PATH:-/tmp/spack-local}
   - export DEVEL_SPACK_CHECKOUT=${DEVEL_SPACK_CHECKOUT:-/tmp/spack}

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Fixed (Repair bugs, etc)
+- [[PR267]](https://github.com/lanl/singularity-eos/pull/267) Add missing eosSafeDestroy call in EOSPAC::Finalize
 - [[PR263]](https://github.com/lanl/singularity-eos/pull/263) Fix stellar collapse mass fraction interpolation
 - [[PR258]](https://github.com/lanl/singularity-eos/pull/258) Fix EOSPAC vector performance and add warnings when slower version gets selected.
 - [[PR243]](https://github.com/lanl/singularity-eos/pull/243) Remove undefined behavior caused by diagnostic variables. Also fixed some compiler warnings.
@@ -13,6 +14,7 @@
 - [[PR232]](https://github.com/lanl/singularity-eos/pull/228) Fixed uninitialized cmake path variables
 
 ### Added (new features/APIs/variables/...)
+- [[PR274]](https://github.com/lanl/singularity-eos/pull/274) added a stiffened gas 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

singularity-eos/eos/eos.hpp

@@ -40,6 +40,7 @@
 #include <singularity-eos/eos/eos_jwl.hpp>
 #include <singularity-eos/eos/eos_spiner.hpp>
 #include <singularity-eos/eos/eos_stellar_collapse.hpp>
+#include <singularity-eos/eos/eos_stiff.hpp>
 #include <singularity-eos/eos/eos_vinet.hpp>
 
 // Modifiers
@@ -63,7 +64,7 @@ using singularity::detail::transform_variadic_list;
 
 // all eos's
 static constexpr const auto full_eos_list =
-    tl<IdealGas, Gruneisen, Vinet, JWL, DavisReactants, DavisProducts
+    tl<IdealGas, Gruneisen, Vinet, JWL, DavisReactants, DavisProducts, StiffGas
 #ifdef SPINER_USE_HDF
        ,
        SpinerEOSDependsRhoT, SpinerEOSDependsRhoSie, StellarCollapse

singularity-eos/eos/eos_eospac.hpp

@@ -1132,7 +1132,10 @@ class EOSPAC : public EosBase<EOSPAC> {
 
   static constexpr unsigned long PreferredInput() { return _preferred_input; }
   PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return 0; }
-  inline void Finalize() {}
+  inline void Finalize() {
+    using namespace EospacWrapper;
+    eosSafeDestroy(NT, tablehandle, Verbosity::Quiet);
+  }
   static std::string EosType() { return std::string("EOSPAC"); }
   static std::string EosPyType() { return EosType(); }
   inline void PrintParams() const { printf("EOSPAC parameters:\nmatid = %s\n", matid_); }

singularity-eos/eos/eos_stiff.hpp

@@ -0,0 +1,201 @@
+//------------------------------------------------------------------------------
+// © 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_STIFF_HPP_
+#define _SINGULARITY_EOS_EOS_EOS_STIFF_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 StiffGas : public EosBase<StiffGas> {
+ public:
+  StiffGas() = default;
+  PORTABLE_INLINE_FUNCTION StiffGas(Real gm1, Real Cv, Real Pinf, Real qq)
+      : _Cv(Cv), _gm1(gm1), _Pinf(Pinf), _qq(qq), _T0(ROOM_TEMPERATURE),
+        _P0(ATMOSPHERIC_PRESSURE), _qp(0.0),
+        _rho0(robust::ratio(_P0 + Pinf, gm1 * Cv * _T0)),
+        _vol0(robust::ratio(gm1 * Cv * _T0, _P0 + Pinf)),
+        _sie0(robust::ratio(_P0 + (gm1 + 1.0) * Pinf, _P0 + Pinf) * Cv * _T0 + qq),
+        _bmod0(gm1 * (gm1 + 1.0) * robust::ratio(_P0 + Pinf, gm1 * Cv * _T0) * Cv * _T0),
+        _dpde0(_rho0 * _gm1) {
+    checkParams();
+  }
+  PORTABLE_INLINE_FUNCTION StiffGas(Real gm1, Real Cv, Real Pinf, Real qq, Real qp,
+                                    Real T0, Real P0)
+      : _Cv(Cv), _gm1(gm1), _Pinf(Pinf), _qq(qq), _T0(T0), _P0(P0), _qp(qp),
+        _rho0(robust::ratio(P0 + Pinf, gm1 * Cv * T0)),
+        _vol0(robust::ratio(gm1 * Cv * _T0, _P0 + Pinf)),
+        _sie0(robust::ratio(P0 + (gm1 + 1.0) * Pinf, P0 + Pinf) * Cv * T0 + qq),
+        _bmod0(gm1 * (gm1 + 1.0) * robust::ratio(P0 + Pinf, gm1 * Cv * T0) * Cv * T0),
+        _dpde0(_rho0 * _gm1) {
+    checkParams();
+  }
+  StiffGas GetOnDevice() { return *this; }
+  PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    return std::max(robust::SMALL(), robust::ratio(rho * (sie - _qq) - _Pinf, rho * _Cv));
+  }
+  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_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    return std::max(_qq, robust::ratio(rho * _Cv * temperature + _Pinf, rho) + _qq);
+  }
+  PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    return std::max(-_Pinf, _gm1 * rho * _Cv * temperature - _Pinf);
+  }
+  PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    return std::max(-_Pinf, _gm1 * rho * (sie - _qq) - (_gm1 + 1.0) * _Pinf);
+  }
+  PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    return _Cv * std::log(robust::ratio(temperature, _T0) + robust::SMALL()) +
+           _gm1 * _Cv * std::log(robust::ratio(_rho0, rho) + robust::SMALL()) + _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 - _Pinf * vol),
+                                        (_sie0 - _qq - _Pinf * _vol0)) +
+                          robust::SMALL()) +
+           _Cv * _gm1 * std::log(robust::ratio(vol, _vol0) + robust::SMALL()) + _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(), _gm1 * (_gm1 + 1.0) * rho * _Cv * temperature);
+  }
+  PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    return std::max(robust::SMALL(), _gm1 * (_gm1 + 1.0) * (rho * (sie - _qq) - _Pinf));
+  }
+  PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature(
+      const Real rho, const Real temperature, Real *lambda = nullptr) const {
+    return _gm1;
+  }
+  PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy(
+      const Real rho, const Real sie, Real *lambda = nullptr) const {
+    return _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(StiffGas)
+  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("Stiff Gas Parameters:\nGamma = %g\nCv    = %g\nPinf  = %g\nq     = "
+           "%g\n",
+           _gm1 + 1.0, _Cv, _Pinf, _qq);
+  }
+  PORTABLE_INLINE_FUNCTION void
+  DensityEnergyFromPressureTemperature(const Real press, const Real temp, Real *lambda,
+                                       Real &rho, Real &sie) const {
+    sie = std::max(
+        _qq,
+        robust::ratio(press + (_gm1 + 1.0) * _Pinf, press + _Pinf) * _Cv * temp + _qq);
+    rho = std::max(robust::SMALL(), robust::ratio(press + _Pinf, _gm1 * _Cv * temp));
+  }
+  inline void Finalize() {}
+  static std::string EosType() { return std::string("StiffGas"); }
+  static std::string EosPyType() { return EosType(); }
+
+ private:
+  Real _Cv, _gm1, _Pinf, _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 StiffGas::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 + (_gm1 + 1.0) * _Pinf, _gm1 * rho) + _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_STIFF_HPP_

singularity-eos/eos/singularity_eos.cpp

@@ -156,6 +156,25 @@ int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0,
                                 Cv0, def_en, def_v);
 }
 
+int init_sg_StiffGas(const int matindex, EOS *eos, const double gm1, const double Cv,
+                     const double Pinf, const double qq, int const *const enabled,
+                     double *const vals) {
+  assert(matindex >= 0);
+  EOS eosi = SGAPPLYMODSIMPLE(StiffGas(gm1, Cv, Pinf, 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(StiffGas(gm1, Cv, Pinf, qq));
+  eos[matindex] = eos_.GetOnDevice();
+  return 0;
+}
+
+int init_sg_StiffGas(const int matindex, EOS *eos, const double gm1, const double Cv,
+                     const double Pinf, const double qq) {
+  return init_sg_StiffGas(matindex, eos, gm1, Cv, Pinf, qq, def_en, def_v);
+}
+
 #ifdef SPINER_USE_HDF
 int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename,
                               const int matid, int const *const enabled,

singularity-eos/eos/singularity_eos.f90

@@ -37,6 +37,7 @@ module singularity_eos
     init_sg_JWL_f,&
     init_sg_DavisProducts_f,&
     init_sg_DavisReactants_f,&
+    init_sg_StiffGas_f,&
     init_sg_SpinerDependsRhoT_f,&
     init_sg_SpinerDependsRhoSie_f,&
     init_sg_eospac_f,&
@@ -120,7 +121,20 @@ end function init_sg_DavisProducts
       type(c_ptr), value, intent(in)         :: sg_mods_enabled, sg_mods_values
     end function init_sg_DavisReactants
   end interface
-
+
+  interface
+    integer(kind=c_int) function &
+      init_sg_StiffGas(matindex, eos, gm1, Cv, Pinf, qq, sg_mods_enabled, &
+                       sg_mods_values) &
+      bind(C, name='init_sg_StiffGas')
+      import
+      integer(c_int), value, intent(in)      :: matindex
+      type(c_ptr), value, intent(in)         :: eos
+      real(kind=c_double), value, intent(in) :: gm1, Cv, Pinf, qq
+      type(c_ptr), value, intent(in)         :: sg_mods_enabled, sg_mods_values
+    end function init_sg_StiffGas
+  end interface
+
   interface
     integer(kind=c_int) function &
       init_sg_SpinerDependsRhoT(matindex, eos, filename, id, sg_mods_enabled, &
@@ -343,7 +357,20 @@ integer function init_sg_DavisReactants_f(matindex, eos, rho0, e0, P0, T0, &
                                  C, G0, Z, alpha, Cv0, c_loc(sg_mods_enabled), &
                                  c_loc(sg_mods_values))
   end function init_sg_DavisReactants_f
-
+
+  integer function init_sg_StiffGas_f(matindex, eos, gm1, Cv, &
+                                      Pinf, qq, &
+                                      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)   :: gm1, Cv, Pinf, qq
+    integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled
+    real(kind=8), dimension(:), target, intent(inout)        :: sg_mods_values
+    err = init_sg_StiffGas(matindex-1, eos%ptr, gm1, Cv, Pinf, qq, &
+                           c_loc(sg_mods_enabled), c_loc(sg_mods_values))
+  end function init_sg_StiffGas_f
+
   integer function init_sg_SpinerDependsRhoT_f(matindex, eos, filename, id, &
                                                sg_mods_enabled, &
                                                sg_mods_values) &

singularity-eos/eos/singularity_eos.hpp

@@ -50,6 +50,10 @@ int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0,
                            const double Cv0, int const *const enabled,
                            double *const vals);
 
+int init_sg_StiffGas(const int matindex, EOS *eos, const double gm1, const double Cv,
+                     const double Pinf, const double qq, int const *const enabled,
+                     double *const vals);
+
 #ifdef SPINER_USE_HDF
 int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename,
                               const int id, int const *const enabled, double *const vals);
@@ -111,6 +115,9 @@ int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0,
                            const double G0, const double Z, const double alpha,
                            const double Cv0);
 
+int init_sg_StiffGas(const int matindex, EOS *eos, const double gm1, const double Cv,
+                     const double Pinf, const double qq);
+
 #ifdef SPINER_USE_HDF
 int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename,
                               const int id);

spack-repo/packages/spiner/package.py

@@ -54,6 +54,7 @@ class Spiner(CMakePackage, CudaPackage):
     depends_on("[email protected]:2.13.9")
     depends_on("[email protected]:", when="@:1.5.1")
     depends_on("[email protected]:", when="@1.6.0:")
+    depends_on("ports-of-call@main", when="@main")
 
     # Currently the raw cuda backend of ports-of-call is not supported.
     depends_on("ports-of-call portability_strategy=Kokkos", when="@:1.5.1 +kokkos")

test/CMakeLists.txt

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