Remove E0 from Davis Products EOS

CHANGELOG.md

@@ -20,10 +20,10 @@
 - [[PR335]](https://github.com/lanl/singularity-eos/pull/335) Fix missing hermite.hpp in CMake install required for Helmholtz EOS
 - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Guard against FPEs in the PTE solver
 - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Update CMake for proper Kokkos linking in Fortran interface
-- [[PR372]](https://github.com/lanl/singularity-eos/pull/372) Fix missing utilization of E0 in Davis Products EOS
 
 ### Changed (changing behavior/API/variables/...)
 - [[PR363]](https://github.com/lanl/singularity-eos/pull/363) Template lambda values for scalar calls
+- [[PR372]](https://github.com/lanl/singularity-eos/pull/372) Removed E0 from Davis Products EOS in favor of using the shifted EOS modifier. CHANGES API!
 
 ### Infrastructure (changes irrelevant to downstream codes)
 - [[PR329]](https://github.com/lanl/singularity-eos/pull/329) Move vinet tests into analytic test suite

doc/sphinx/src/models.rst

@@ -1231,7 +1231,7 @@ Finally, the pressure, energy, and temperature along the isentrope are given by
 
 .. math::
 
-    e_S(\rho) = e_{\mathrm{C}} G(\rho) \frac{1}{\rho V_{\mathrm{C}}} - e_0
+    e_S(\rho) = e_{\mathrm{C}} G(\rho) \frac{1}{\rho V_{\mathrm{C}}}
 
 .. math::
 
@@ -1272,25 +1272,29 @@ The energy at the CJ state can be calculated as
     e_\mathrm{CJ} = \frac{P_0 + P_\mathrm{CJ}}{2(V_0 - V_\mathrm{CJ})},
 
 relative to :math:`e = 0` at the reference state of the *reactants*. Therefore
-the :math:`e_0` energy offset of the products EOS is given by
+the energy offset of the products EOS is given by
 
 .. math::
 
     e_0 = e_S(V_\mathrm{CJ}) - e_\mathrm{CJ}.
 
 Practically, this means :math:`e_0` should be positive for any energetic material.
 
+To provide the energy offset to the Davis Products EOS, `the energy shift
+modifier<modifiers shifted EOS>`_ should be used. Note that the convention there
+is that the shift is positive, so :math:`-e_0` should be provided to the shift
+modifier.
+
 The constructor for the Davis Products EOS is
 
 .. code-block:: cpp
 
   DavisProducts(const Real a, const Real b, const Real k, const Real n, const Real vc,
-                const Real pc, const Real Cv, const Real E0)
+                const Real pc, const Real Cv)
 
 where ``a`` is :math:`a`, ``b`` is :math:`b`, ``k`` is :math:`k`,
 ``n`` is :math:`n`, ``vc`` is :math:`V_\mathrm{C}`, ``pc`` is
-:math:`P_\mathrm{C}`, ``Cv`` is :math:`C_{V,0}`, and ``E0`` is
-:math:`e_\mathrm{C}`.
+:math:`P_\mathrm{C}`, ``Cv`` is :math:`C_{V,0}`.
 
 Spiner EOS
 ````````````

doc/sphinx/src/modifiers.rst

@@ -22,6 +22,8 @@ We list below the available modifiers and their constructors.
 The Shifted EOS
 -----------------
 
+.. _modifiers shifted EOS:
+
 The shifted equation of state modifies zero point energy of an
 underlying model by some shift. So for example, it transforms
 

python/module.cpp

@@ -62,7 +62,7 @@ PYBIND11_MODULE(singularity_eos, m) {
     .def(
       py::init<Real, Real, Real, Real, Real, Real, Real, Real>(),
       py::arg("a"), py::arg("b"), py::arg("k"), py::arg("n"), py::arg("vc"),
-      py::arg("pc"), py::arg("Cv"), py::arg("E0")
+      py::arg("pc"), py::arg("Cv")
     );
 
 #ifdef SPINER_USE_HDF

singularity-eos/eos/eos_davis.hpp

@@ -177,8 +177,8 @@ class DavisProducts : public EosBase<DavisProducts> {
   DavisProducts() = default;
   PORTABLE_INLINE_FUNCTION
   DavisProducts(const Real a, const Real b, const Real k, const Real n, const Real vc,
-                const Real pc, const Real Cv, const Real E0)
-      : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv), _E0(E0) {}
+                const Real pc, const Real Cv)
+      : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv) {}
   DavisProducts GetOnDevice() { return *this; }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy(
@@ -286,8 +286,8 @@ class DavisProducts : public EosBase<DavisProducts> {
   static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; }
   PORTABLE_INLINE_FUNCTION void PrintParams() const {
     static constexpr char s1[]{"DavisProducts Params: "};
-    printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e E0:%e\n", s1, _a, _b, _k, _n, _vc,
-           _pc, _Cv, _E0);
+    printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e \n", s1, _a, _b, _k, _n, _vc, _pc,
+           _Cv);
   }
   inline void Finalize() {}
   static std::string EosType() { return std::string("DavisProducts"); }
@@ -296,7 +296,7 @@ class DavisProducts : public EosBase<DavisProducts> {
 
  private:
   static constexpr Real onethird = 1.0 / 3.0;
-  Real _a, _b, _k, _n, _vc, _pc, _Cv, _E0;
+  Real _a, _b, _k, _n, _vc, _pc, _Cv;
   // static constexpr const char _eos_type[] = "DavisProducts";
   static constexpr const unsigned long _preferred_input =
       thermalqs::density | thermalqs::specific_internal_energy;
@@ -314,8 +314,7 @@ class DavisProducts : public EosBase<DavisProducts> {
     const Real ec = _pc * _vc / (_k - 1.0 + _a);
     // const Real de = ecj-(Es(rho0)-_E0);
     return ec * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) /
-               std::pow(vvc, _k - 1.0 + _a) -
-           _E0;
+           std::pow(vvc, _k - 1.0 + _a);
   }
   PORTABLE_INLINE_FUNCTION Real Ts(const Real rho) const {
     const Real vvc = 1 / (rho * _vc);

singularity-eos/eos/singularity_eos.cpp

@@ -154,23 +154,23 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B,
 
 int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b,
                           const double k, const double n, const double vc,
-                          const double pc, const double Cv, const double E0,
-                          int const *const enabled, double *const vals) {
+                          const double pc, const double Cv, int const *const enabled,
+                          double *const vals) {
   assert(matindex >= 0);
-  EOS eosi = SGAPPLYMODSIMPLE(DavisProducts(a, b, k, n, vc, pc, Cv, E0));
+  EOS eosi = SGAPPLYMODSIMPLE(DavisProducts(a, b, k, n, vc, pc, Cv));
   if (enabled[3] == 1) {
     singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2],
                                            vals[3], vals[4], vals[5]);
   }
-  EOS eos_ = SGAPPLYMOD(DavisProducts(a, b, k, n, vc, pc, Cv, E0));
+  EOS eos_ = SGAPPLYMOD(DavisProducts(a, b, k, n, vc, pc, Cv));
   eos[matindex] = eos_.GetOnDevice();
   return 0;
 }
 
 int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b,
                           const double k, const double n, const double vc,
-                          const double pc, const double Cv, const double E0) {
-  return init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, E0, def_en, def_v);
+                          const double pc, const double Cv) {
+  return init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, def_en, def_v);
 }
 
 int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0,

singularity-eos/eos/singularity_eos.f90

@@ -112,13 +112,13 @@ end function init_sg_JWL
 
   interface
     integer(kind=c_int) function &
-      init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, E0, &
+      init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, &
                             sg_mods_enabled, sg_mods_values) &
       bind(C, name='init_sg_DavisProducts')
       import
       integer(c_int), value, intent(in)      :: matindex
       type(c_ptr), value, intent(in)         :: eos
-      real(kind=c_double), value, intent(in) :: a, b, k, n, vc, pc, Cv, E0
+      real(kind=c_double), value, intent(in) :: a, b, k, n, vc, pc, Cv
       type(c_ptr), value, intent(in)         :: sg_mods_enabled, sg_mods_values
     end function init_sg_DavisProducts
   end interface
@@ -460,12 +460,12 @@ integer function init_sg_JWL_f(matindex, eos, A, B, R1, R2, w, rho0, Cv, &
   end function init_sg_JWL_f
 
   integer function init_sg_DavisProducts_f(matindex, eos, a, b, k, n, vc, pc, &
-                                           Cv, E0, sg_mods_enabled, &
+                                           Cv, 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)   :: a, b, k, n, vc, pc, Cv, E0
+    real(kind=8), value, intent(in)   :: a, b, k, n, vc, pc, Cv
     integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled
     real(kind=8), dimension(:), target, optional, intent(inout)        :: sg_mods_values
     ! local vars
@@ -478,7 +478,7 @@ integer function init_sg_DavisProducts_f(matindex, eos, a, b, k, n, vc, pc, &
     if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values
 
     err = init_sg_DavisProducts(matindex-1, eos%ptr, a, b, k, n, vc, pc, Cv, &
-                                E0, c_loc(sg_mods_enabled_use), &
+                                c_loc(sg_mods_enabled_use), &
                                 c_loc(sg_mods_values_use))
   end function init_sg_DavisProducts_f
 

singularity-eos/eos/singularity_eos.hpp

@@ -39,8 +39,8 @@ int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const doubl
 
 int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b,
                           const double k, const double n, const double vc,
-                          const double pc, const double Cv, const double E0,
-                          int const *const enabled, double *const vals);
+                          const double pc, const double Cv, int const *const enabled,
+                          double *const vals);
 
 int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0,
                            const double e0, const double P0, const double T0,
@@ -134,7 +134,7 @@ int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const doubl
 
 int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b,
                           const double k, const double n, const double vc,
-                          const double pc, const double Cv, const double E0);
+                          const double pc, const double Cv);
 
 int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0,
                            const double e0, const double P0, const double T0,

test/pte_test_utils.hpp

@@ -59,7 +59,7 @@ inline void set_eos(T *eos) {
   singularity::EOS dr = singularity::DavisReactants(
       1.890, 4.115e10, 1.0e6, 297.0, 1.8e5, 4.6, 0.34, 0.56, 0.0, 0.4265, 0.001074e10);
   singularity::EOS dp = singularity::DavisProducts(0.798311, 0.58, 1.35, 2.66182, 0.75419,
-                                                   3.2e10, 0.001072e10, 0.0);
+                                                   3.2e10, 0.001072e10);
   eos[0] = gr.GetOnDevice();
   eos[1] = dr.GetOnDevice();
   eos[2] = dp.GetOnDevice();

test/test_f_iface.f90

@@ -43,7 +43,7 @@ program test_sg_fortran_interface
                             4.6d0, 0.34d0, 0.56d0,0.d0, 0.4265d0, 0.001074d10)
 mat = mat + 1
 res = init_sg_DavisProducts_f(mat, eos, 0.798311d0, 0.58d0, 1.35d0, 2.66182d0,&
-                           0.75419d0, 3.2d10, 0.001072d10, 0.d0)
+                           0.75419d0, 3.2d10, 0.001072d10)
 
 ! cleanup
 res = finalize_sg_eos_f(nmat, eos)