Skip to content

Commit

Permalink
Merge pull request #372 from lanl/jhp/DavisEOS
Browse files Browse the repository at this point in the history
Remove E0 from Davis Products EOS
  • Loading branch information
jhp-lanl authored May 14, 2024
2 parents 994da55 + 768f629 commit a213659
Show file tree
Hide file tree
Showing 11 changed files with 57 additions and 27 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
33 changes: 30 additions & 3 deletions doc/sphinx/src/models.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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<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
4 changes: 2 additions & 2 deletions python/module.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -60,9 +60,9 @@ PYBIND11_MODULE(singularity_eos, m) {
eos_class<DavisProducts>(m, "DavisProducts")
.def(py::init())
.def(
py::init<Real, Real, Real, Real, Real, Real, Real, Real>(),
py::init<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
10 changes: 5 additions & 5 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,16 +286,16 @@ 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"); }
static std::string EosPyType() { return EosType(); }

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 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
6 changes: 3 additions & 3 deletions singularity-eos/eos/singularity_eos.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand Down
2 changes: 1 addition & 1 deletion test/pte_test_utils.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down
2 changes: 1 addition & 1 deletion test/python_bindings.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
2 changes: 1 addition & 1 deletion test/test_f_iface.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit a213659

Please sign in to comment.