diff --git a/CHANGELOG.md b/CHANGELOG.md index 0220ebcb58..465399f4cf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,7 @@ ### 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 diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 12f555caac..ca42caed77 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1254,17 +1254,44 @@ Here, there are four dimensionless parameters that are settable by the user, :math:`e_\mathrm{C}`, :math:`V_\mathrm{C}` and :math:`T_\mathrm{C}` are tuning parameters with units related to their non-subscripted counterparts. +Note that the energy zero (i.e. the reference energy) for the Davis products EOS +is arbitrary. For the isentrope to properly pass through the CJ state of a +reacting material, the energy release of the reaction needs to be accounted for +properly. If done external to the EOS, an energy source term is required in the +Euler equations. However, a common convention is to specify the reactants and +product EOS in a consistent way such that the reference energy corresponds to +the rest state of the material *before* it reacts. + +The energy at the CJ state can be calculated as + +.. math:: + + 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 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`_ 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 ```````````` diff --git a/doc/sphinx/src/modifiers.rst b/doc/sphinx/src/modifiers.rst index 552b16940b..75e11310af 100644 --- a/doc/sphinx/src/modifiers.rst +++ b/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 diff --git a/python/module.cpp b/python/module.cpp index a6b3ad1c79..2b89795552 100644 --- a/python/module.cpp +++ b/python/module.cpp @@ -60,9 +60,9 @@ PYBIND11_MODULE(singularity_eos, m) { eos_class(m, "DavisProducts") .def(py::init()) .def( - py::init(), + py::init(), 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 diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 87e2ef06f6..0fb2217e0d 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -177,8 +177,8 @@ class DavisProducts : public EosBase { 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 PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( @@ -286,8 +286,8 @@ class DavisProducts : public EosBase { 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"); } @@ -295,7 +295,7 @@ class DavisProducts : public EosBase { 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; diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index e000bfe268..c80ab710b5 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/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, diff --git a/singularity-eos/eos/singularity_eos.f90 b/singularity-eos/eos/singularity_eos.f90 index 2c0ef227ec..e325601758 100644 --- a/singularity-eos/eos/singularity_eos.f90 +++ b/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 diff --git a/singularity-eos/eos/singularity_eos.hpp b/singularity-eos/eos/singularity_eos.hpp index 9a7ddd43f7..9a31d10f54 100644 --- a/singularity-eos/eos/singularity_eos.hpp +++ b/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, diff --git a/test/pte_test_utils.hpp b/test/pte_test_utils.hpp index d61b09747e..af2e68ad94 100644 --- a/test/pte_test_utils.hpp +++ b/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(); diff --git a/test/python_bindings.py b/test/python_bindings.py index a0ff479ad7..951e90523e 100644 --- a/test/python_bindings.py +++ b/test/python_bindings.py @@ -87,7 +87,7 @@ def testDavisReactants(self): eos = singularity_eos.DavisReactants(1,1,1,1,1,1,1,1,1,1,1) def testDavisProducts(self): - eos = singularity_eos.DavisProducts(1,1,1,1,1,1,1,1) + eos = singularity_eos.DavisProducts(1,1,1,1,1,1,1) class Modifiers(unittest.TestCase, EOSTestBase): def setUp(self): diff --git a/test/test_f_iface.f90 b/test/test_f_iface.f90 index 52a92e95cc..7b845a2119 100644 --- a/test/test_f_iface.f90 +++ b/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)