Vinet EOS implemented

doc/sphinx/src/models.rst

@@ -326,6 +326,78 @@ There is an overload of the ``Gruneisen`` class which computes
   Gruneisen(const Real C0, const Real s1, const Real s2, const Real s3, const Real G0,
             const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv)
 
+Extender Vinet EOS
+``````````````````
+
+The extended Vinet EOS is a full EOS, extended in both temperature and density from the Vinet 
+universal EOS for solids (also called Rose cold curve). It is expected to work well in 
+compression but is untested in expansion. It is published in Appendix 2 
+in `J. Appl. Phys. 119, 015904 (2016)`_.
+
+.. _J. Appl. Phys. 119, 015904 (2016): https://doi.org/10.1063/1.4939675
+
+While the Mie-Gruneisen EOS is based on a Hugoniot as reference curve, the Vinet is based on an isotherm:
+
+.. math::
+
+    P(\rho,T) = P_{ref}(\rho) + \alpha_0 B_0 (T - T_{ref})
+
+where the reference isotherm is
+
+.. math::
+
+    P_{ref}(X)&=\frac{3 B_0}{X^2} Z \exp[\eta_0 Z] \left( 1 + \sum_{n=2}^N d_n Z^n \right) \, , \\
+    X &= \left( \frac{\rho_0}{\rho} \right)^{1/3} \\
+    Z &= 1-X
+
+Note that :math:`P_{ref}=0` when :math:`\rho = \rho_0`, the reference state on the reference
+isotherm is always at ambient pressure. However, the reference isotherm is not necessarily
+at room temperature.
+
+It can be shown that :math:`B_0` is the isothermal bulk modulus, and :math:`\alpha_0` the 
+thermal expansion coefficient, at the reference state, and that
+
+.. math::
+
+    \eta_0 = \frac{3}{2}\left[ \left[ \frac{\partial B}{\partial P}\right]_0 -1\right] \, .
+
+By assuming that also the constant volume heat capacity is a constant, :math:`{C_V}_0`, an entropy
+can be derived
+
+.. math:: 
+
+    S(V,T) = S_0 + \alpha_0 B_0 (V - V_0) + {C_V}_0 \ln \frac{T}{T_{ref}}
+
+and from that a thermodynamic consistent energy
+
+.. math::
+
+  E(X,T) =& 9 \frac{B_0 V_0}{{\eta_0}^2}\left(f_0 - \exp[\eta_0 Z] \left(f_0 - \eta_0 Z \left(f_0 + \sum_{n=1}^N f_n Z^n \right)\right)\right) \\
+          & - \alpha_0 B_0 V_0 (1-X^3) T_{ref} + {C_V}_0 (T - T_{ref}) + E_0
+
+where the energy coefficients :math:`f_n` are determined from the pressure coefficients :math:`d_n`, :math:`n\geq 2`, by
+
+.. math::
+
+    f_N &= d_N \\
+    f_n &= d_n - \frac{n+2}{\eta_0} f_{n+1} \\
+    d_0 &= 1.0 \\
+    d_1 &= 0.0 
+
+
+The constructor for the ``Vinet`` EOS has the signature
+
+.. code-block:: cpp
+
+ Vinet(const Real rho0, const Real T0, const Real B0, const Real BP0, const Real A0,
+            const Real Cv0, const Real E0, const Real S0, const Real *expconsts)
+
+where ``rho0`` is :math:`\rho_0`, ``T0`` is :math:`T_{ref}`, ``B0`` is
+:math:`B_0`, ``BP0`` is :math:`(\partial B/\partial P)_0`, ``A0`` is :math:`\alpha_0`, 
+``Cv0`` is :math:`{C_V}_0`, ``E0`` is :math:`E_0`, ``S0`` is :math:`S_0`, and
+``expconsts`` is a pointer to the constant array of length 39 containing the expansion coefficients
+:math:`d_2` to :math:`d_{40}`. Expansion coefficients not used should be set to 0.0.
+
 JWL EOS
 ````````
 

singularity-eos/CMakeLists.txt

@@ -24,6 +24,7 @@ set(EOS_HEADERS
     eos/eos_davis.hpp
     eos/eos_eospac.hpp
     eos/eos_gruneisen.hpp
+    eos/eos_vinet.hpp
     eos/eos_ideal.hpp
     eos/eos_jwl.hpp
     eos/eos_spiner.hpp

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_vinet.hpp>
 
 // Modifiers
 #include <singularity-eos/eos/modifiers/eos_unitsystem.hpp>
@@ -62,7 +63,7 @@ using singularity::detail::transform_variadic_list;
 
 // all eos's
 static constexpr const auto full_eos_list =
-    tl<IdealGas, Gruneisen, JWL, DavisReactants, DavisProducts
+    tl<IdealGas, Gruneisen, Vinet, JWL, DavisReactants, DavisProducts
 #ifdef SPINER_USE_HDF
        ,
        SpinerEOSDependsRhoT, SpinerEOSDependsRhoSie, StellarCollapse

singularity-eos/eos/eos_builder.cpp

@@ -196,6 +196,18 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par
     Real Cv = mpark::get<Real>(base_params["Cv"]);
     return Gruneisen(C0, s1, s2, s3, G0, b, rho0, T0, P0, Cv);
   }
+  //  if (type == EOSType::Vinet) {
+  //    Real rho0 = mpark::get<Real>(base_params["rho0"]);
+  //    Real T0 = mpark::get<Real>(base_params["T0"]);
+  //    Real B0 = mpark::get<Real>(base_params["B0"]);
+  //    Real BP0 = mpark::get<Real>(base_params["BP0"]);
+  //    Real A0 = mpark::get<Real>(base_params["A0"]);
+  //    Real Cv0 = mpark::get<Real>(base_params["Cv0"]);
+  //    Real E0 = mpark::get<Real>(base_params["E0"]);
+  //    Real S0 = mpark::get<Real>(base_params["S0"]);
+  //    Real d2to40 = mpark::get<Real>(base_params["d2to40"]);
+  //    return Vinet(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40);
+  //  }
   if (type == EOSType::JWL) {
     GETBASE(Real, A); // I got tired of writing this over and over
     GETBASE(Real, B);

singularity-eos/eos/eos_builder.hpp

@@ -33,6 +33,7 @@ namespace EOSBuilder {
 enum class EOSType {
   IdealGas,
   Gruneisen,
+  //  Vinet,
   JWL,
   DavisReactants,
   DavisProducts