Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Remove E0 from Davis Products EOS #372

Merged
merged 18 commits into from
May 14, 2024
Merged
Show file tree
Hide file tree
Changes from 10 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 9 additions & 5 deletions doc/sphinx/src/models.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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::

Expand Down Expand Up @@ -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
````````````
Expand Down
2 changes: 2 additions & 0 deletions doc/sphinx/src/modifiers.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
2 changes: 1 addition & 1 deletion python/module.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
13 changes: 6 additions & 7 deletions singularity-eos/eos/eos_davis.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -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"); }
Expand All @@ -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;
Expand All @@ -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);
Expand Down
12 changes: 6 additions & 6 deletions singularity-eos/eos/singularity_eos.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
10 changes: 5 additions & 5 deletions singularity-eos/eos/singularity_eos.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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

Expand Down
Loading