diff --git a/.gitignore b/.gitignore index 08aaff4dcf..d4dc8e6dc9 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,6 @@ utils/variant cmake-gen* spack-build-* spack-configure-* +._* +*.sublime* + diff --git a/CHANGELOG.md b/CHANGELOG.md index ff1c433dca..3b0a01c643 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,12 +2,11 @@ ## Current develop -### Fixed (Repair, bugs, ect) +### Fixed (Repair bugs, etc) - [[PR238]](https://github.com/lanl/singularity-eos/pull/238) Fixed broken examples -- [[PR232]](https://github.com/lanl/singularity-eos/pull/228) Fixed uninitialized cmake path variables -- [[PR228]](https://github.com/lanl/singularity-eos/pull/228) and [[PR229]](https://github.com/lanl/singularity-eos/pull/229) added untracked header files in cmake - [[PR228]](https://github.com/lanl/singularity-eos/pull/228) added untracked header files in cmake - [[PR215]](https://github.com/lanl/singularity-eos/pull/215) and [[PR216]](https://github.com/lanl/singularity-eos/pull/216) fix duplicate definition of EPS and fix CI +- [[PR232]](https://github.com/lanl/singularity-eos/pull/228) Fixed uninitialized cmake path variables ### Added (new features/APIs/variables/...) - [[PR226]](https://github.com/lanl/singularity-eos/pull/226) added entropy interpolation to stellar collapse eos @@ -15,12 +14,15 @@ - [[PR226]](https://github.com/lanl/singularity-eos/pull/226) added entropy interpolation to stellar collapse eos - [[PR209]](https://github.com/lanl/singularity-eos/pull/209) added more documentation around how to contribute to the project and also what a contributor can expect from the core team - [[PR214]](https://github.com/lanl/singularity-eos/pull/214) added documentation about adding a new EOS +- [[PR228]](https://github.com/lanl/singularity-eos/pull/228) and [[PR229]](https://github.com/lanl/singularity-eos/pull/229) added untracked header files in cmake +- [[PR233]](https://github.com/lanl/singularity-eos/pull/233) Added entropy infrastructure to singularity. Add entropy for the ideal gas and modifiers while throwing an error for EOS where entropy is not implemented yet - [[PR177]](https://github.com/lanl/singularity-eos/pull/177) added EOSPAC vector functions ### Changed (changing behavior/API/variables/...) - [[PR223]](https://github.com/lanl/singularity-eos/pull/223) Update ports-of-call and add portable error handling - [[PR234]](https://github.com/lanl/singularity-eos/pull/234) update ports-of-call to correct for undefined behavior in error handling - [[PR219]](https://github.com/lanl/singularity-eos/pull/219) Removed static analysis from re-git pipeline +- [[PR233]](https://github.com/lanl/singularity-eos/pull/233) Exposed entropy for the EOS type (now required for future EOS) ### Infrastructure (changes irrelevant to downstream codes) - [[PR190]](https://github.com/lanl/singularity-eos/pull/190) update CI on re-git diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 0084113b7d..f9b15d3c0f 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1,5 +1,21 @@ .. _models: +.. + Important references: + +.. _MenikoffAndPlohr: https://doi.org/10.1103/RevModPhys.61.75 + +.. _MaxwellWiki: https://en.wikipedia.org/wiki/Maxwell_relations + +.. _MenikoffCompleteEOS: https://www.osti.gov/biblio/1241653 + +.. _DavisReactants: https://doi.org/10.1016/S0010-2180(99)00112-1 + +.. _ProbingOffHugoniotStates: https://doi.org/10.1063/1.4939675 + +.. _WillsThermo: https://www.osti.gov/biblio/1561015 + + EOS Models =========== @@ -7,7 +23,7 @@ The mathematical descriptions of these models are presented below while the details of using them is presented in the description of the :doc:`EOS API `. -EOS Primer +EOS Theory ---------- An equation of state (EOS) is a constituitive model that generally relates @@ -22,15 +38,13 @@ information needed to derive a given thermodynamic quantity. However, not all variables are exposed since the emphasis is on only those that are needed for hydrodynamics codes. An incomplete equation of state is often sufficient for pure hydrodynamics since it can relate pressure to a given energy-density state, -but it is missing information that relates the energy to temperature (which in -turn provides a means to calculate entropy). +but it is missing thermal information (see the section on :ref:`complete +equations of state `). An excellent resource on equations of state in the context of hydrodynamics is -the seminal work from `Menikoff and Plohr`_. In particular, Appendix A contains -a number of thermodynamic relationships that can be useful for computing -additional quantities from those output in ``singularity-eos``. - -.. _Menikoff and Plohr: https://doi.org/10.1103/RevModPhys.61.75 +the seminal work from `Menikoff and Plohr `_. In particular, +Appendix A contains a number of thermodynamic relationships that can be useful +for computing additional quantities from those output in ``singularity-eos``. The Mie-Gruneisen form ```````````````````````` @@ -59,51 +73,206 @@ Davis EOS use the material's isentrope. In ths way, the reference curve indicates the conditions under which you can expect the EOS to represent the intended behavior. -Nomenclature ---------------------- +Some Notes on Thermodynamic Consistency +```````````````````````````````````````` + +.. _Complete EOS: + +For the pure purpose of solving the Euler equations, an incomplete equation of +state of the form :math:`P(\rho, e)` is sufficient. In essence, this provides the +mechanical response of a material subjected to different types of compression +and expansion. However, this EOS is lacking *thermal* information without an +appropriate heat capacity relationship. + +As discussed `by Mattsson `_, an equation of state can be +considered "complete" if it is constructed from one of the thermodynamic +potentials using their *natural variables*, i.e. + +.. math:: + + e =& e(\rho, S) \\ + h =& h(P, S) \\ + F =& F(\rho, T) \\ + G =& G(P, T) \\ + +where all potentials are specific. Here :math:`e` is again the internal energy, +:math:`h` is the enthalpy, :math:`F` is the Helmholtz free energy, and :math:`G` +is the Gibbs free energy. While equations of state formulated using the +Helmholtz free energy can be particularly attractive (such as the sesame +tables), finding a convenient form can be difficult. As such, it becomes +imperitive to extend the Mie-Gruneisen form so that it can form a complete +EOS. + +The heat capacity is defined as + +.. math:: + + C_V := \left(\frac{\partial e}{\partial T}\right)_V + = T \left(\frac{\partial S}{\partial T}\right)_V + +which provides a natural means by which to relate the entropy/energy to the +temperature and form a complete equation of state. However, there are +specific requirements placed on the functional form of the heat capacity +described below. + + +The `Maxwell relations `_ are the consequence of the +requirement that mixed partial second derivatives of the thermodynamic +potentials should be independent of the order of differentiation. This concept +can be extended to third derivatives (as Davis does in his Appendix B `here +`_) to produce the relationship, + +.. math:: + + \frac{V}{C_V^2}\left(\frac{\partial C_V}{\partial V}\right)_S = + \left(\frac{\partial \Gamma}{\partial S}\right)_V. + +This is often referred to as a "compatibility condition" (see also +`Menikoff `_) and provides an important connection +between the Gruneisen parameter, :math:`\Gamma`, and the constant volume heat +capacity, :math:`C_V`, in the context of a complete equation of state +developed from the internal energy. Importantly, if the Gruneisen form forces +us to assume that the Gruneisen parameter is a sole function of +density/volume, then the implication is that the heat capacity ***must*** then +be a sole function of entropy. + +Again from `Menikoff `_, it can be shown that even +in the simplest possible case of a constant heat capacity, the temperature +is related to the energy through + +.. math:: + + T(\rho, e) = T_0 \phi(\rho) + \frac{e - e_{S,0}(\rho)}{C_V} + +where :math:`T_0` represents a reference temperature, :math:`e_{S,0}` is the +energy along the isentrope that passes through the reference temperature, and +:math:`\phi(\rho)` is an integrating factor given by + +.. math:: + + \phi(\rho) = \exp\left(\int\limits_{\rho_0}^{\rho} \rho \Gamma + (\rho) \mathrm{d}\rho \right). + +As an EOS of a Mie-Gruneisen form becomes more complicated with more complex +functional forms for :math:`\Gamma` and the reference curves, the task of +calculating a **thermodynamically consistent** temperature becomes more +complicated. + +Available EOS Information and Nomenclature +------------------------------------------ The EOS models in ``singularity-eos`` are defined for the following sets of dependent and independent variables through various member functions described in the :doc:`EOS API `. -+--------------------------+----------------------+--------------------------+ -| Function | Dependent Variable | Independent Variables | -+==========================+======================+==========================+ -| :math:`T(\rho, e)` | Temperature | Density, Internal Energy | -+--------------------------+----------------------+ | -| :math:`P(\rho, e)` | Pressure | | -+--------------------------+----------------------+--------------------------+ -| :math:`e(\rho, T)` | Internal Energy | Density, Temperature | -+--------------------------+----------------------+ | -| :math:`P(\rho, T)` | Pressure | | -+--------------------------+----------------------+--------------------------+ -| :math:`\rho(P, T)` | Density | Pressure, Temperature | -+--------------------------+----------------------+ | -| :math:`e(P, T)` | Internal Energy | | -+--------------------------+----------------------+--------------------------+ -| :math:`C_V(\rho, T)` | Constant Volume | Density, Temperature | -+--------------------------+ Specific Heat +--------------------------+ -| :math:`C_V(\rho, e)` | Capacity | Density, Internal Energy | -+--------------------------+----------------------+--------------------------+ -| :math:`B_S(\rho, T)` | Isentropic Bulk | Density, Temperature | -+--------------------------+ Modulus +--------------------------+ -| :math:`B_S(\rho, e)` | | Density, Internal Energy | -+--------------------------+----------------------+--------------------------+ -| :math:`\Gamma(\rho, T)` | Gruneisen Parameter | Density, Temperature | -+--------------------------+ +--------------------------+ -| :math:`\Gamma(\rho, e)` | | Density, Internal Energy | -+--------------------------+----------------------+--------------------------+ ++---------------------------+----------------------+--------------------------+ +| Function | Dependent Variable | Independent Variables | ++===========================+======================+==========================+ +| :math:`T(\rho, e)` | Temperature | Density, Specific | ++---------------------------+----------------------+ Internal Energy | +| :math:`P(\rho, e)` | Pressure | | ++---------------------------+----------------------+ | +| :math:`S(\rho, e)` | Specific Entropy | | ++---------------------------+----------------------+--------------------------+ +| :math:`e(\rho, T)` | Specific Internal | Density, Temperature | +| | Energy | | ++---------------------------+----------------------+ | +| :math:`P(\rho, T)` | Pressure | | ++---------------------------+----------------------+ | +| :math:`S(\rho, T)` | Specific Entropy | | ++---------------------------+----------------------+--------------------------+ +| :math:`\rho(P, T)` [#PT]_ | Density | Pressure, [#PT]_ | ++---------------------------+----------------------+ Temperature [#PT]_ | +| :math:`e(P, T)` [#PT]_ | Specific Internal | | +| | Energy | | ++---------------------------+----------------------+--------------------------+ +| :math:`C_V(\rho, T)` | Constant Volume | Density, Temperature | ++---------------------------+ Specific Heat +--------------------------+ +| :math:`C_V(\rho, e)` | Capacity | Density, Specific | +| | | Internal Energy | ++---------------------------+----------------------+--------------------------+ +| :math:`B_S(\rho, T)` | Isentropic Bulk | Density, Temperature | ++---------------------------+ Modulus +--------------------------+ +| :math:`B_S(\rho, e)` | | Density, Specific | +| | | Internal Energy | ++---------------------------+----------------------+--------------------------+ +| :math:`\Gamma(\rho, T)` | Gruneisen Parameter | Density, Temperature | ++---------------------------+ +--------------------------+ +| :math:`\Gamma(\rho, e)` | | Density, Specific | +| | | Internal Energy | ++---------------------------+----------------------+--------------------------+ + +.. [#PT] + Note: Using pressure and temperature as independent variables is fraught + since both pressure and energy are often multi-valued in density for many + EOS due to the presence of phase changes (especially tabular EOS). For + analytic EOS in ``singularity-eos`` where there is not an analytic inversion + to pressure-temperature space, a root-find is typically used that uses the + density at standard temperature and pressure (STP) as an initial guess. The + notable exceptions to this are the spiner EOS that allow the introduction of + an initial guess via the ``lambda`` function argument. A point of note is that "specific" implies that the quantity is intensive on a -per unit mass basis. It should be assumed that the internal energy is *always* -specific since we are working in terms of density (the inverse of specific -volume). +per unit mass basis. It should be assumed that the internal energy and entopry +are *always* specific since we are working in terms of density (the inverse of +specific volume). -Disambiguation -```````````````` +Entropy availability +```````````````````` +For an arbitrary equation of state, a change in entropy in terms of temperature +and volume is given by -Gruneisen Parameter -''''''''''''''''''' +.. math:: + + \Delta S = \int\limits_{S_0}^{S} \mathrm{d}S = + \int\limits_{T_0}^{T} \left(\frac{\partial S}{\partial T}\right)_V + \mathrm{d}T + + \int\limits_{V_0}^{V} \left(\frac{\partial S}{\partial V}\right)_T + \mathrm{d}V, + +which can be simplified using a definition of the heat capacity and a Maxwell +relation to become + +.. math:: + + \Delta S = + \int\limits_{T_0}^{T} \frac{C_V}{T} \mathrm{d}T + + \int\limits_{V_0}^{V} \left(\frac{\partial P}{\partial T}\right)_V + \mathrm{d}V. + +Similarly, expressing the entropy in terms of *energy* and volume yields + +.. math:: + + \Delta S = + \int\limits_{e_0}^{e} \frac{1}{T(\rho, e)} \mathrm{d}e + + \int\limits_{V_0}^{V} \frac{P(\rho, e)}{T(\rho, e)} + \mathrm{d}V + +after substituting for the appropriate derivatives using the first law of +thermodynamics. + +Importantly for the analytic EOS, these integrals will tend to diverge as the +temperature and volume approach zero if the heat capacity does not also +approach zero. This necessitates appropriate choices for the reference +states :math:`T_0` and :math:`V_0`. + +.. note:: + + All EOS objects will expose the functions + :code:`EntropyFromDensityInternalEnergy ()` and + :code:`EntropyFromDensityTemperature()` but many EOS cannot currently + calculate entropy, either because the EOS is not thermodynamically + consistent or because the feature has not yet been implemented. In these + cases, the use of these functions will abort the calculation or raise an + exception depending on the execution space (host or device). The cases where + this occurs are noted below. + +Nomenclature Disambiguation +```````````````````````````` + +The Gruneisen Parameter +''''''''''''''''''''''' In this description of the EOS models, we use :math:`\Gamma` to represent the Gruneisen coeficient since this is the most commonly-used symbol in the context of Mie-Gruneisen equations of state. The definition of the Gruneisen @@ -147,6 +316,8 @@ units for thermodynamic quantities: +--------------+------------------+---------------------------------------+-----------------------+ |:math:`e` | erg/g | 10\ :sup:`-12` Terg/g | 10\ :sup:`-10` J/mg | +--------------+------------------+---------------------------------------+-----------------------+ +|:math:`S` | erg/g-K | 10\ :sup:`-12` Terg/g-K | 10\ :sup:`-10` J/mg-K | ++--------------+------------------+---------------------------------------+-----------------------+ |:math:`T` | K | 1 | 1 | +--------------+------------------+---------------------------------------+-----------------------+ |:math:`C_V` | erg/g/K | 10\ :sup:`-12` Terg/g/K | 10\ :sup:`-10` J/mg/K | @@ -156,7 +327,7 @@ units for thermodynamic quantities: |:math:`\Gamma`| unitless | -- | -- | +--------------+------------------+---------------------------------------+-----------------------+ -Note: some temperatures are measured in eV for which the conversion is +Note: sometimes temperatures are measured in eV for which the conversion is 8.617333262e-05 eV/K. Sesame units are equivalent to the mm-mg-µs unit system. @@ -179,25 +350,57 @@ the form e = C_V T, -where quantities are defined in the :ref:`nomenclature ` section. +where quantities are defined in the :ref:`nomenclature ` section. A common point of confusion +is :math:`\Gamma` versus :math:`\gamma` with the latter being the adiabatic +exponent. For an ideal gas, they are related through -The settable parameters are the Gruneisen parameter and specific heat capacity. +.. math:: + \Gamma = \gamma - 1 -Although this differs from the traditional representation of the ideal gas law -as :math:`P\tilde{V} = RT`, the forms are equivalent by recognizing that -:math:`\Gamma = \frac{R}{\tilde{C_V}}` where :math:`R` is the ideal gas constant -in units of energy per mole per Kelvin and :math:`\tilde{C_\mathrm{V}}` is the -*molar* heat capacity, relatable to the *specific* heat capacity through the -molecular weight of the gas. Since :math:`\tilde{C_\mathrm{V}} = \frac{5}{2} R` -for a diatomic ideal gas, the corresponding Gruneisen parameter should be 0.4. +Although this formulation differs from the traditional representation of the +ideal gas law as :math:`P\tilde{V} = RT`, the forms are equivalent by +recognizing that :math:`\Gamma = \frac{R}{\tilde{C_V}}` where :math:`R` is the +ideal gas constant in units of energy per mole per Kelvin and :math:`\tilde +{C_\mathrm{V}}` is the *molar* heat capacity, relatable to the *specific* heat +capacity through the molecular weight of the gas. Since :math:`\tilde{C_\mathrm +{V}} = \frac{5}{2} R` for a diatomic ideal gas, the corresponding Gruneisen +parameter should be 0.4. -A common point of confusion is :math:`\Gamma` versus :math:`\gamma` with the -latter being the adiabatic exponent. For an ideal gas, they are related through +The entropy for an ideal gas is given by .. math:: - \Gamma = \gamma - 1 + S = C_V \ln\left(\frac{T}{T_0}\right) + \Gamma C_V \ln\left(\frac{\rho_0} + {\rho}\right), + +.. note:: + + The entropy diverges to negative infinity at absolute zero due to the + constant heat capacity assumption. Care should be taken when using + temperatures significantly below that of the reference state. + +we have assumed that the entropy is zero at the reference state given +by :math:`T_0` and :math:`\rho_0`. By default, :math:`T_0 = 298` K and the +reference density is given by + +.. math:: + + \rho_0 = \frac{P_0}{\Gamma C_V T_0}, + +where :math:`P_0 = 1` bar. + +It should be noted that this equation diverges as the temperature approaches +zero or the density approaches infinity. From a physical perspective, this is a +limitation of the constant heat capacity assumption and would be remedied if +the heat capacity went to zero at absolute zero. However, this represents a +significant deviation from true ideal gas physics so we do not include it +here. + +The settable parameters for this EOS are the Gruneisen parameter and specific +heat capacity. Optionally, the reference state for the entropy calculation can +be provided by setting *both* the reference density and temperature. The ``IdealGas`` EOS constructor has two arguments, ``gm1``, which is the Gruneisen parameter :math:`\Gamma`, and ``Cv``, which is the @@ -205,11 +408,27 @@ specific heat :math:`C_V`: .. code-block:: cpp - IdealGas(Real gm1, Real Cv) + IdealGas(Real gm1, Real Cv) + +Optionally, the reference temperature and density for the entropy calculation, +can be provided in the constructor via ``EntropyT0`` and ``EntropyRho0``: + +.. code-block:: cpp + + IdealGas(Real gm1, Real Cv, Real EntropyT0, Real EntropyRho0) + +Note that these parameters are provided solely for the entropy calculation. When +these values are not set, they will be the same as those returned by the +:code:`ValuesAtReferenceState()` function. However, if the entropy reference +conditions are given, the return values of the :code:`ValuesAtReferenceState()` +function will not be the same. Gruneisen EOS ````````````` +.. warning:: + Entropy is not available for this EOS + One of the most commonly-used EOS to represent solids is the Steinberg variation of the Mie-Gruneisen EOS, often just shortened to "Gruneisen" EOS. This EOS uses the Hugoniot as the reference curve and thus is extremly powerful because @@ -222,8 +441,8 @@ The pressure follows the traditional Mie-Gruneisen form, P(\rho, e) = P_H(\rho) + \rho\Gamma(\rho) \left(e - e_H(\rho) \right), Here the subscript :math:`H` is a reminder that the reference curve is a -Hugoniot. Other quantities are defined in the :ref:`nomenclature ` -section. +Hugoniot. Other quantities are defined in the :ref:`nomenclature ` section. The above is an incomplete equation of state because it only relates the pressure to the density and energy, the minimum required in a solution to the @@ -238,6 +457,10 @@ The user should note that this implies that :math:`e=0` at the reference temperature, :math:`T_0`. Given this simple relationship, the user should treat the temperature from this EOS as only a rough estimate. +Given the inconsisetency in the temperature, we have made the choice **not** to +expose the entropy for this EOS. **Requesting an entropy value will result in an +error.** + Given a reference density, :math:`\rho_0`, we first parameterize the EOS using :math:`\eta` as a measure of compression given by @@ -285,7 +508,8 @@ where :math:`P_0` is the reference pressure and :math:`c_0`, :math:`s_1`, Here :math:`U_s` is the shock velocity and :math:`u_p` is the particle velocity. For many materials, this relationship is roughly linear so only the :math:`s_1` parameter is needed. The units for :math:`c_0` are velocity while -the rest are unitless. +the rest are unitless. Note that the parameter :math:`s_1` is related to the +fundamental derivative of shock physics as shown by `Wills `_. Finally the energy along the Hugoniot is given by @@ -326,17 +550,17 @@ 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 -`````````````````` +Extendended 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 +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) `_. -While the Mie-Gruneisen EOS is based on a Hugoniot as reference curve, the Vinet is based on an isotherm: +While the Mie-Gruneisen EOS is based on a Hugoniot as reference curve, the Vinet +is based on an isotherm: .. math:: @@ -350,19 +574,20 @@ where the reference isotherm is 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. +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 +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 +By assuming that also the constant volume heat capacity is a constant, +:math:`{C_V}_0`, an entropy can be derived .. math:: @@ -372,18 +597,25 @@ 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) \\ + 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 +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 + d_1 &= 0.0 + +.. note:: + The entropy diverges to negative infinity at absolute zero due to the + constant heat capacity assumption. Care should be taken when using + temperatures significantly below that of the reference state. The constructor for the ``Vinet`` EOS has the signature @@ -393,13 +625,18 @@ The constructor for the ``Vinet`` EOS has the signature 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. +: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 -```````` +`````````` + +.. warning:: + Entropy is not available for this EOS The Jones-Wilkins-Lee (JWL) EOS is used mainly for detonation products of high explosives. Similar to the other EOS here, the JWL EOS can be written in a @@ -446,12 +683,13 @@ where ``A`` is :math:`A`, ``B`` is :math:`B`, ``R1`` is :math:`R_1`, and ``Cv`` is :math:`C_V`. Davis EOS -````````` +``````````` -The Davis reactants and products EOS are both of Mie-Gruneisen forms that use -isentropes for the reference curves. The equations of state are typically used -to represent high explosives and their detonation products and the reference -curves are calibrated to several sets of experimental data. +The `Davis reactants `_ and products EOS are both of +Mie-Gruneisen forms that use isentropes for the reference curves. The equations +of state are typically used to represent high explosives and their detonation +products and the reference curves are calibrated to several sets of +experimental data. For both the reactants and products EOS, the pressure and energy take the forms @@ -465,51 +703,97 @@ For both the reactants and products EOS, the pressure and energy take the forms \right), where the subscript :math:`S` denotes quantities along the reference isentrope -and other quantities are defined in the :ref:`nomenclature ` -section. +and other quantities are defined in the :ref:`nomenclature ` section. Davis Reactants EOS ''''''''''''''''''' -The Davis reactants EOS uses an isentrope passing through a reference state -and assumes that the heat capacity varies linearly with entropy such that +.. warning:: + Entropy is not yet available for this EOS + +The `Davis reactants EOS `_ uses an isentrope passing through a +reference state and as the reference curve and then assumes that the heat +capacity varies linearly with entropy such that .. math:: C_V = C_{V,0} + \alpha(S - S_0), where subscript :math:`0` refers to the reference state and :math:`\alpha` is -a dimensionless constant specified by the user. +a dimensionless constant specified by the user. -The :math:`e(\rho, P)` lookup is quite awkward, so the energy is -more-conveniently cast in terms of termperature such that +The Gruneisen parameter is given a linear form such that .. math:: - e(\rho, T) = e_S(\rho) + \frac{C_{V,0} T_S(\rho)}{1 + \alpha} - \left( \left(\frac{T}{T_S(\rho)} \right)^{1 + \alpha} - 1 \right), + \Gamma(\rho) = \Gamma_0 + + \begin{cases} + 0 & \rho < \rho_0 \\ + Zy & \rho >= \rho_0 + \end{cases} -which can easily be inverted to find :math:`T(\rho, e)`. +where :math:`Z` is a dimensionless parameter and :math:`y = 1 - \rho0/\rho`. +Along an isentrope, the Gruneisen parameter can be expressed as + +.. math:: + + \Gamma_S(\rho) = \frac{\rho}{T} + \left(\frac{\partial T}{\partial \rho}\right)_S, -The Gruneisen parameter takes on a linear form such that +which, upon integration can produce the temperature along the reference +isentrope: .. math:: - \Gamma(\rho) = \Gamma_0 + + T_{S,0}(\rho) = \begin{cases} - 0 & \rho < \rho_0 \\ - Zy & \rho >= \rho_0 + T_0\left(\frac{\rho}{\rho_0}\right)^{\Gamma_0} & \rho < \rho_0 \\ + T_0\exp\left(-Zy\right)\left(\frac{\rho}{\rho_0}\right)^{\Gamma_0 + Z} + & \rho >= \rho_0 \end{cases} -where :math:`Z` and :math:`y` are dimensionless parameters. +where :math:`T_{S,0}` is the temperature along the reference isentrope, +:math:`S = S_0`. -Finally, the pressure, energy, and temperature along the isentrope are given by +Using the fact that the heat capacity can be expressed as + +.. math:: + + C_V = T\left( \frac{\partial S}{\partial T} \right)_V, + +the temperature off of the reference isoentrope can be integrated from this +identity to yield + +.. math:: + + T(\rho, S) = T_{S,0}(\rho) \left( \frac{C_V(S)}{C_{V,0}} \right)^{\frac{1}{\alpha}}, + +Now requiring that the entropy go to zero at absolute zero in accordance with +the Nernst postulate and the third law of thermodynamics, the entropy can be +expressed as a function of temperature and density such that + +.. math:: + + S(\rho, T) = \frac{C_{V,0}}{\alpha} \left( \frac{T}{T_{S,0}(\rho)} \right)^\alpha. + +The :math:`e(\rho, P)` formulation can now be more-conveniently cast in terms of +termperature such that + +.. math:: + + e(\rho, T) = e_S(\rho) + \frac{C_{V,0} T_S(\rho)}{1 + \alpha} + \left( \left(\frac{T}{T_S(\rho)} \right)^{1 + \alpha} - 1 \right), + +which can easily be inverted to find :math:`T(\rho, e)`. + +Finally, the pressure and energy along the isentrope are given by .. math:: P_S(\rho) = P_0 + \frac{\rho_0 A^2}{4B} \begin{cases} - \mathrm{e} \left( 4By \right) -1 & \rho < \rho_0 \\ + \exp \left( 4By \right) -1 & \rho < \rho_0 \\ \sum\limits_{j=1}^3 \frac{(4By)^j}{j!} + C\frac{(4By)^4}{4!} + \frac{y^2}{(1-y)^4} & \rho >= \rho0 \end{cases} @@ -517,16 +801,8 @@ Finally, the pressure, energy, and temperature along the isentrope are given by .. math:: e_S(\rho) = e_0 + \int\limits_{\rho_0}^{\rho} - \frac{P_S(\bar{\rho})}{\bar{\rho^2}}~\mathrm{d}\bar{\rho} - -.. math:: + \frac{P_S(\bar{\rho})}{\bar{\rho^2}}~\mathrm{d}\bar{\rho}, - T_S(\rho) = T_0 - \begin{cases} - \left(\frac{\rho}{\rho_0} \right)^{\Gamma_0} & \rho < \rho_0 \\ - \mathrm{e} \left( -Zy \right) \left(\frac{\rho}{\rho_0} \right)^{\Gamma_0 + Z} - & \rho >= \rho_0 - \end{cases} where :math:`A`, :math:`B`, :math:`C`, :math:`y`, and :math:`Z` are all user-settable parameters and again quantities with a subcript of :math:`0` @@ -556,6 +832,9 @@ heat capacity at the reference state. Davis Products EOS ''''''''''''''''''' +.. warning:: + Entropy is not yet available for this EOS + The Davis products EOS is created from the reference isentrope passing through the CJ state of the high explosive along with a constant heat capacity. The constant heat capacity leads to the energy being a simple funciton of the @@ -631,6 +910,9 @@ where ``a`` is :math:`a`, ``b`` is :math:`b`, ``k`` is :math:`k`, Spiner EOS ```````````` +.. warning:: + Entropy is not yet available for this EOS + Spiner EOS is a tabulated reader for the `Sesame`_ database of material equations of state. Materials include things like water, dry air, iron, or steel. This model comes in two flavors: @@ -842,6 +1124,9 @@ which saves the current EOS data in ``sp5`` format. EOSPAC EOS ```````````` +.. warning:: + Entropy is not yet available for this EOS + This is a striaghtforward wrapper of the `EOSPAC`_ library for the `Sesame`_ database. The constructor for the ``EOSPAC`` model looks like diff --git a/doc/sphinx/src/modifiers.rst b/doc/sphinx/src/modifiers.rst index 5060557340..5eb2e78151 100644 --- a/doc/sphinx/src/modifiers.rst +++ b/doc/sphinx/src/modifiers.rst @@ -65,7 +65,44 @@ energy. The pressure is unchanged under the operation for some scale parameter :math:`s`. The ``ScaledEOS`` applies this transformation to any equation of state, not just an ideal gas, where -the pressure may change. +the pressure may change for different scaling ratios. + +Another way of understanding scaling ratios is that the pressure can be written +as + +.. math:: + + P = \left(\frac{\partial F}{\partial V} \right)_T + +where :math:`F` is the Helmholtz free energy. For a given scaling such that +:math:`\rho_\mathrm{eos} = s\rho_\mathrm{in}`, the volume obeys the inverse +scaling. Since the scaling ratio is constant, it can be substituted into the +above expression so that + +.. math:: + + P = \left(\frac{\partial F_\mathrm{eos}}{\partial V_\mathrm{eos}} \right)_T + = \left(\frac{\partial F_\mathrm{in}}{\partial V_\mathrm{in}} \right)_T + = \left(\frac{\partial F_\mathrm{in}}{s \partial V_\mathrm{eos}} \right)_T + = \left(\frac{s\partial F_\mathrm{eos}}{s \partial V_\mathrm{eos}} \right)_T + +which implies that the Helmholtz free energy must scale in the same way as +volume (inverse to density) in order to preserve the same pressure. Applying +this scaling to the definition of the Helmholtz free energy yields + +.. math:: + + F_\mathrm{eos} = e_\mathrm{eos} - TS_\mathrm{eos} = \frac{1}{R} F_\mathrm{in} + = \frac{1}{R}e_\mathrm{in} - T\left(\frac{1}{R}S_\mathrm{in}\right), + +where the implicaiton is that this inverse the scaling ratio should also be +applied to energy. The inverse scaling ratio must be applied to the entropy +here in order to ensure that all other thermodynamic potentials +(energy, entropy, and the Gibbs free energy) scale similarly. + +where :math:`e` is the internal energy and :math:`S` is the entropy. The +implication is that the same scaling should be applied to the energy and entropy +to maintain thermodynamic consistency. The constructor for ``ScaledEOS`` accepts the underlying model, and the scale parameter. For example, a shifted ideal gas might be diff --git a/doc/sphinx/src/sphinx-doc.rst b/doc/sphinx/src/sphinx-doc.rst index 4764209519..d0c347139b 100644 --- a/doc/sphinx/src/sphinx-doc.rst +++ b/doc/sphinx/src/sphinx-doc.rst @@ -33,15 +33,14 @@ with shared machines or python distributions. First, ensure that you are running a modern version of python (i.e. python 3 of some flavor) -.. code-block:: - :language:bash +.. code-block:: bash + $ python --version Python 3.9.7 Then, use pip to install :code:`spinx` and the RTD theme -.. code-block:: - :language:bash +.. code-block:: bash pip install --user sphinx sphinx-rtd-theme @@ -49,6 +48,7 @@ Now, navigate to the :code:`../doc/sphinx` directory where a :code:`make help` shows all of the available ways to build the documentation .. code-block:: + :language:bash $ make help diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index f481d51e0c..26b3faaef9 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -15,9 +15,11 @@ #ifndef _SINGULARITY_EOS_EOS_EOS_BASE_ #define _SINGULARITY_EOS_EOS_EOS_BASE_ +#include #include #include +#include namespace singularity { namespace mfuncname { @@ -50,7 +52,10 @@ namespace eos_base { using EosBase::MinimumDensity; \ using EosBase::MinimumTemperature; \ using EosBase::PTofRE; \ - using EosBase::FillEos; + using EosBase::FillEos; \ + using EosBase::EntropyFromDensityTemperature; \ + using EosBase::EntropyFromDensityInternalEnergy; \ + using EosBase::EntropyIsNotEnabled; /* This is a CRTP that allows for static inheritance so that default behavior for @@ -150,6 +155,48 @@ class EosBase { PressureFromDensityInternalEnergy(rhos, sies, pressures, num, lambdas); } template + inline void EntropyFromDensityTemperature(ConstRealIndexer &&rhos, + ConstRealIndexer &&temperatures, + RealIndexer &&entropies, const int num, + LambdaIndexer &&lambdas) const { + static auto const name = SG_MEMBER_FUNC_NAME(); + static auto const cname = name.c_str(); + CRTP copy = *(static_cast(this)); + portableFor( + cname, 0, num, PORTABLE_LAMBDA(const int i) { + entropies[i] = + copy.EntropyFromDensityTemperature(rhos[i], temperatures[i], lambdas[i]); + }); + } + template + inline void + EntropyFromDensityTemperature(ConstRealIndexer &&rhos, ConstRealIndexer &&temperatures, + RealIndexer &&entropies, Real * /*scratch*/, + const int num, LambdaIndexer &&lambdas) const { + EntropyFromDensityTemperature(rhos, temperatures, entropies, num, lambdas); + } + template + inline void EntropyFromDensityInternalEnergy(ConstRealIndexer &&rhos, + ConstRealIndexer &&sies, + RealIndexer &&entropies, const int num, + LambdaIndexer &&lambdas) const { + static auto const name = SG_MEMBER_FUNC_NAME(); + static auto const cname = name.c_str(); + CRTP copy = *(static_cast(this)); + portableFor( + cname, 0, num, PORTABLE_LAMBDA(const int i) { + entropies[i] = + copy.EntropyFromDensityInternalEnergy(rhos[i], sies[i], lambdas[i]); + }); + } + template + inline void + EntropyFromDensityInternalEnergy(ConstRealIndexer &&rhos, ConstRealIndexer &&sies, + RealIndexer &&entropies, Real * /*scratch*/, + const int num, LambdaIndexer &&lambdas) const { + EntropyFromDensityInternalEnergy(rhos, sies, entropies, num, lambdas); + } + template inline void SpecificHeatFromDensityTemperature(ConstRealIndexer &&rhos, ConstRealIndexer &&temperatures, RealIndexer &&cvs, const int num, @@ -339,69 +386,15 @@ class EosBase { PORTABLE_INLINE_FUNCTION Real RhoPmin(const Real temp) const { return 0.0; } - // Specialzied vector version of PTofRE maybe more suited for EOSPAC - // } - // template - // inline - // void PTofRE(RealIndexer &&rhos, RealIndexer &&sies, - // RealIndexer &&presses, RealIndexer &&temps, - // RealIndexer &&dpdrs, RealIndexer &&dpdes, - // RealIndexer &&dtdrs, RealIndexer &&dtdes, - // const int num, LambdaIndexer &&lambdas) const { - // // Get the dervived class - // auto eos = static_cast(*this); - - // // Lookup at density and energy - // eos.PressureFromDensityInternalEnergy(rhos, sies, num, presses, - // lambdas); - // eos.TemperatureFromDensityInternalEnergy(rhos, sies, num, temps, - // lambdas); - // // Peturbation factor - // Real factor = 1. + 1.0e-06; - - // // Perturb densities and do lookups - // portableFor( - // 'PerturbDensities', 0, num, PORTABLE_LAMBDA(const int i) { - // rhos[i] *= factor; - // } - // ); - // eos.PressureFromDensityInternalEnergy(rhos, sies, num, dpdrs, - // lambdas); - // eos.TemperatureFromDensityInternalEnergy(rhos, sies, num, dtdrs, - // lambdas); - - // // Reset densities, perturb energies, and do lookups - // portableFor( - // 'PerturbEnergiesResetDensities', 0, num, - // PORTABLE_LAMBDA(const int i) { - // sies[i] *= factor; - // rhos[i] /= factor; - // } - // ); - // eos.PressureFromDensityInternalEnergy(rhos, sies, dpdes, num, - // lambdas); - // eos.TemperatureFromDensityInternalEnergy(rhos, sies, dtdes, num, - // lambdas); - - // // Reset the energies to their original values - // portableFor( - // 'ResetEnergies', 0, num, PORTABLE_LAMBDA(const int i) { - // sies[i] /= factor; - // } - // ); - - // // Calculate the derivatives - // portableFor( - // 'CalculateDerivatives', 0., num, PORTABLE_LAMBDA(const int i) { - // dpdrs[i] = (dpdrs[i] - presses[i]) / (rhos[i] * (1. - factor)); - // dpdes[i] = (dpdes[i] - presses[i]) / (sies[i] * (1. - factor)); - // dtdrs[i] = (dtdrs[i] - temps[i]) / (rhos[i] * (1. - factor)); - // dtdes[i] = (dtdes[i] - temps[i]) / (sies[i] * (1. - factor)); - // } - // ); - - // return; - // } + // Default entropy behavior is to return an error + PORTABLE_FORCEINLINE_FUNCTION + void EntropyIsNotEnabled(const char *eosname) const { + // Construct the error message using char* so it works on device + char msg[120] = "Entropy is not enabled for the '"; + std::strcat(msg, eosname); + std::strcat(msg, "' EOS"); + PORTABLE_ALWAYS_THROW_OR_ABORT(msg); + } }; } // namespace eos_base } // namespace singularity diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 37d7c31a43..c6eaa1f3f3 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -60,6 +60,16 @@ class DavisReactants : public EosBase { const Real rho, const Real sie, Real *lambda = nullptr) const { return Ps(rho) + Gamma(rho) * rho * (sie - Es(rho)); } + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + EntropyIsNotEnabled("DavisReactants"); + return 1.0; + } + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + EntropyIsNotEnabled("DavisReactants"); + return 1.0; + } PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( const Real rho, const Real temp, Real *lambda = nullptr) const { return SpecificHeatFromDensityInternalEnergy( @@ -152,6 +162,16 @@ class DavisProducts : public EosBase { const Real rho, const Real sie, Real *lambda = nullptr) const { return Ps(rho) + rho * Gamma(rho) * (sie - Es(rho)); } + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + EntropyIsNotEnabled("DavisProducts"); + return 1.0; + } + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + EntropyIsNotEnabled("DavisProducts"); + return 1.0; + } PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const { return _Cv; diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp index 2d5f13657e..721e3c93a1 100644 --- a/singularity-eos/eos/eos_eospac.hpp +++ b/singularity-eos/eos/eos_eospac.hpp @@ -54,6 +54,10 @@ class EOSPAC : public EosBase { const Real rho, const Real temperature, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( @@ -82,6 +86,8 @@ class EOSPAC : public EosBase { using EosBase::InternalEnergyFromDensityTemperature; using EosBase::PressureFromDensityTemperature; using EosBase::PressureFromDensityInternalEnergy; + using EosBase::EntropyFromDensityTemperature; + using EosBase::EntropyFromDensityInternalEnergy; using EosBase::SpecificHeatFromDensityTemperature; using EosBase::SpecificHeatFromDensityInternalEnergy; using EosBase::BulkModulusFromDensityTemperature; @@ -90,6 +96,7 @@ class EOSPAC : public EosBase { using EosBase::GruneisenParamFromDensityInternalEnergy; using EosBase::PTofRE; using EosBase::FillEos; + using EosBase::EntropyIsNotEnabled; // EOSPAC vector implementations template @@ -810,6 +817,12 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::PressureFromDensityTemperature(const Real return Real(pressureFromSesame(P[0])); } +PORTABLE_INLINE_FUNCTION Real EOSPAC::EntropyFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda) const { + EntropyIsNotEnabled("EOSPAC"); + return 1.0; +} + PORTABLE_INLINE_FUNCTION void EOSPAC::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, const unsigned long output, @@ -923,6 +936,12 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::PressureFromDensityInternalEnergy( Real temp = TemperatureFromDensityInternalEnergy(rho, sie, lambda); return PressureFromDensityTemperature(rho, temp, lambda); } +PORTABLE_INLINE_FUNCTION Real EOSPAC::EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda) const { + using namespace EospacWrapper; + const Real temp = TemperatureFromDensityInternalEnergy(rho, sie, lambda); + return EntropyFromDensityTemperature(rho, temp, lambda); +} PORTABLE_INLINE_FUNCTION Real EOSPAC::SpecificHeatFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda) const { using namespace EospacWrapper; diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index 237a5f35f9..ece535f0d5 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -76,6 +76,10 @@ class Gruneisen : public EosBase { const Real rho, const Real temperature, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( const Real rho, const Real temperatummmmmmre, Real *lambda = nullptr) const { return _Cv; @@ -310,6 +314,12 @@ PORTABLE_INLINE_FUNCTION Real Gruneisen::PressureFromDensityInternalEnergy( } return P_H + Gamma(rho) * rho * (sie - E_H); } +PORTABLE_INLINE_FUNCTION Real Gruneisen::EntropyFromDensityInternalEnergy( + const Real rho_in, const Real sie, Real *lambda) const { + const Real rho = std::min(rho_in, _rho_max); + EntropyIsNotEnabled("Gruneisen"); + return 1.0; +} PORTABLE_INLINE_FUNCTION Real Gruneisen::BulkModulusFromDensityInternalEnergy( const Real rho_in, const Real sie, Real *lambda) const { using namespace gruneisen_utils; @@ -332,6 +342,12 @@ PORTABLE_INLINE_FUNCTION Real Gruneisen::PressureFromDensityTemperature( return PressureFromDensityInternalEnergy( rho, InternalEnergyFromDensityTemperature(rho, temp)); } +PORTABLE_INLINE_FUNCTION Real Gruneisen::EntropyFromDensityTemperature( + const Real rho_in, const Real temp, Real *lambda) const { + const Real rho = std::min(rho_in, _rho_max); + const Real sie = InternalEnergyFromDensityTemperature(rho, temp); + return EntropyFromDensityInternalEnergy(rho, sie); +} PORTABLE_INLINE_FUNCTION Real Gruneisen::BulkModulusFromDensityTemperature( const Real rho_in, const Real temp, Real *lambda) const { const Real rho = std::min(rho_in, _rho_max); diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp index a281c091e7..27af667d31 100644 --- a/singularity-eos/eos/eos_ideal.hpp +++ b/singularity-eos/eos/eos_ideal.hpp @@ -26,6 +26,7 @@ // Base stuff #include #include +#include #include #define MYMAX(a, b) a > b ? a : b @@ -40,13 +41,31 @@ class IdealGas : public EosBase { PORTABLE_INLINE_FUNCTION IdealGas(Real gm1, Real Cv) : _Cv(Cv), _gm1(gm1), _rho0(_P0 / (_gm1 * _Cv * _T0)), _sie0(_Cv * _T0), _bmod0((_gm1 + 1) * _gm1 * _rho0 * _Cv * _T0), _dpde0(_gm1 * _rho0), - _dvdt0(1. / (_rho0 * _T0)) {} + _dvdt0(1. / (_rho0 * _T0)), _EntropyT0(_T0), _EntropyRho0(_rho0) { + checkParams(); + } + PORTABLE_INLINE_FUNCTION IdealGas(Real gm1, Real Cv, Real EntropyT0, Real EntropyRho0) + : _Cv(Cv), _gm1(gm1), _rho0(_P0 / (_gm1 * _Cv * _T0)), _sie0(_Cv * _T0), + _bmod0((_gm1 + 1) * _gm1 * _rho0 * _Cv * _T0), _dpde0(_gm1 * _rho0), + _dvdt0(1. / (_rho0 * _T0)), _EntropyT0(EntropyT0), _EntropyRho0(EntropyRho0) { + checkParams(); + } IdealGas GetOnDevice() { return *this; } PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda = nullptr) const { return MYMAX(0.0, sie / _Cv); } + PORTABLE_INLINE_FUNCTION void checkParams() const { + // Portable_require seems to do the opposite of what it should. Conditions + // reflect this and the code should be changed when ports-of-call changes + PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Heat capacity must be positive"); + PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive"); + PORTABLE_ALWAYS_REQUIRE(_EntropyT0 >= 0, + "Entropy reference temperature must be positive"); + PORTABLE_ALWAYS_REQUIRE(_EntropyRho0 >= 0, + "Entropy reference density must be positive"); + } PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const { return MYMAX(0.0, _Cv * temperature); @@ -59,6 +78,16 @@ class IdealGas : public EosBase { const Real rho, const Real sie, Real *lambda = nullptr) const { return MYMAX(0.0, _gm1 * rho * sie); } + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + return _Cv * log(robust::ratio(temperature, _EntropyT0)) + + _gm1 * _Cv * log(robust::ratio(_EntropyRho0, rho)); + } + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + const Real temp = TemperatureFromDensityInternalEnergy(rho, sie, lambda); + return EntropyFromDensityTemperature(rho, temp, lambda); + } PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const { return _Cv; @@ -132,6 +161,8 @@ class IdealGas : public EosBase { // static constexpr const char _eos_type[] = {"IdealGas"}; static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; + // optional entropy reference state variables + Real _EntropyT0, _EntropyRho0; }; PORTABLE_INLINE_FUNCTION diff --git a/singularity-eos/eos/eos_jwl.hpp b/singularity-eos/eos/eos_jwl.hpp index 17b14e62ec..fd0ad207c0 100644 --- a/singularity-eos/eos/eos_jwl.hpp +++ b/singularity-eos/eos/eos_jwl.hpp @@ -50,6 +50,10 @@ class JWL : public EosBase { const Real rho, const Real temperature, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( @@ -118,6 +122,12 @@ PORTABLE_INLINE_FUNCTION Real JWL::PressureFromDensityInternalEnergy(const Real Real *lambda) const { return ReferencePressure(rho) + _w * rho * (sie - ReferenceEnergy(rho)); } +PORTABLE_INLINE_FUNCTION Real JWL::EntropyFromDensityInternalEnergy(const Real rho, + const Real sie, + Real *lambda) const { + EntropyIsNotEnabled("JWL"); + return 1.0; +} PORTABLE_INLINE_FUNCTION Real JWL::TemperatureFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda) const { return (sie - ReferenceEnergy(rho)) / _Cv; @@ -146,6 +156,12 @@ PORTABLE_INLINE_FUNCTION Real JWL::PressureFromDensityTemperature(const Real rho return PressureFromDensityInternalEnergy( rho, InternalEnergyFromDensityTemperature(rho, temp)); } +PORTABLE_INLINE_FUNCTION Real JWL::EntropyFromDensityTemperature(const Real rho, + const Real temp, + Real *lambda) const { + EntropyIsNotEnabled("JWL"); + return 1.0; +} PORTABLE_INLINE_FUNCTION Real JWL::SpecificHeatFromDensityTemperature( const Real rho, const Real temp, Real *lambda) const { return SpecificHeatFromDensityInternalEnergy( diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index 60242febc4..9dbd138a84 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -77,6 +77,8 @@ class SpinerEOSDependsRhoT : public EosBase { using EosBase::InternalEnergyFromDensityTemperature; using EosBase::PressureFromDensityTemperature; using EosBase::PressureFromDensityInternalEnergy; + using EosBase::EntropyFromDensityTemperature; + using EosBase::EntropyFromDensityInternalEnergy; using EosBase::SpecificHeatFromDensityTemperature; using EosBase::SpecificHeatFromDensityInternalEnergy; using EosBase::BulkModulusFromDensityTemperature; @@ -85,6 +87,7 @@ class SpinerEOSDependsRhoT : public EosBase { using EosBase::GruneisenParamFromDensityInternalEnergy; using EosBase::PTofRE; using EosBase::FillEos; + using EosBase::EntropyIsNotEnabled; inline SpinerEOSDependsRhoT(const std::string &filename, int matid, bool reproduciblity_mode = false); @@ -109,6 +112,12 @@ class SpinerEOSDependsRhoT : public EosBase { Real PressureFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION + Real EntropyFromDensityTemperature(const Real rho, const Real temperature, + Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION + Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(const Real rho, const Real temperature, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION @@ -300,6 +309,8 @@ class SpinerEOSDependsRhoSie : public EosBase { using EosBase::InternalEnergyFromDensityTemperature; using EosBase::PressureFromDensityTemperature; using EosBase::PressureFromDensityInternalEnergy; + using EosBase::EntropyFromDensityTemperature; + using EosBase::EntropyFromDensityInternalEnergy; using EosBase::SpecificHeatFromDensityTemperature; using EosBase::SpecificHeatFromDensityInternalEnergy; using EosBase::BulkModulusFromDensityTemperature; @@ -308,6 +319,7 @@ class SpinerEOSDependsRhoSie : public EosBase { using EosBase::GruneisenParamFromDensityInternalEnergy; using EosBase::PTofRE; using EosBase::FillEos; + using EosBase::EntropyIsNotEnabled; PORTABLE_INLINE_FUNCTION SpinerEOSDependsRhoSie() : memoryStatus_(DataStatus::Deallocated) {} @@ -333,6 +345,12 @@ class SpinerEOSDependsRhoSie : public EosBase { Real PressureFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION + Real EntropyFromDensityTemperature(const Real rho, const Real temperature, + Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION + Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(const Real rho, const Real T, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION @@ -932,6 +950,21 @@ Real SpinerEOSDependsRhoT::PressureFromDensityInternalEnergy(const Real rho, return P; } +PORTABLE_INLINE_FUNCTION +Real SpinerEOSDependsRhoT::EntropyFromDensityTemperature(const Real rho, + const Real temperature, + Real *lambda) const { + EntropyIsNotEnabled("SpinerEOSDependsRhoT"); + return 1.0; +} +PORTABLE_INLINE_FUNCTION +Real SpinerEOSDependsRhoT::EntropyFromDensityInternalEnergy(const Real rho, + const Real sie, + Real *lambda) const { + EntropyIsNotEnabled("SpinerEOSDependsRhoT"); + return 1.0; +} + PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::SpecificHeatFromDensityTemperature(const Real rho, const Real temperature, @@ -1648,6 +1681,21 @@ Real SpinerEOSDependsRhoSie::PressureFromDensityInternalEnergy(const Real rho, return interpRhoSie_(rho, sie, dependsRhoSie_.P, lambda); } +PORTABLE_INLINE_FUNCTION +Real SpinerEOSDependsRhoSie::EntropyFromDensityTemperature(const Real rho, + const Real temperature, + Real *lambda) const { + EntropyIsNotEnabled("SpinerEOSDependsRhoSie"); + return 1.0; +} +PORTABLE_INLINE_FUNCTION +Real SpinerEOSDependsRhoSie::EntropyFromDensityInternalEnergy(const Real rho, + const Real sie, + Real *lambda) const { + EntropyIsNotEnabled("SpinerEOSDependsRhoSie"); + return 1.0; +} + PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoSie::SpecificHeatFromDensityTemperature(const Real rho, const Real T, diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 338ae252d6..75282ba9ec 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -71,6 +71,8 @@ class StellarCollapse : public EosBase { using EosBase::InternalEnergyFromDensityTemperature; using EosBase::PressureFromDensityTemperature; using EosBase::PressureFromDensityInternalEnergy; + using EosBase::EntropyFromDensityTemperature; + using EosBase::EntropyFromDensityInternalEnergy; using EosBase::SpecificHeatFromDensityTemperature; using EosBase::SpecificHeatFromDensityInternalEnergy; using EosBase::BulkModulusFromDensityTemperature; @@ -107,6 +109,9 @@ class StellarCollapse : public EosBase { Real EntropyFromDensityTemperature(const Real rho, const Real temperature, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION + Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(const Real rho, const Real temperature, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION @@ -503,6 +508,15 @@ Real StellarCollapse::EntropyFromDensityTemperature(const Real rho, return (entropy > robust::EPS() ? entropy : robust::EPS()); } +PORTABLE_INLINE_FUNCTION +Real StellarCollapse::EntropyFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda) const { + Real lRho, lT, Ye; + getLogsFromRhoSie_(rho, sie, lambda, lRho, lT, Ye); + const Real entropy = entropy_.interpToReal(Ye, lT, lRho); + return (entropy > robust::EPS() ? entropy : robust::EPS()); +} + PORTABLE_INLINE_FUNCTION Real StellarCollapse::SpecificHeatFromDensityTemperature(const Real rho, const Real temperature, diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index a2d97772f4..2bd2c86447 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -111,6 +111,25 @@ class Variant { }, eos_); } + PORTABLE_INLINE_FUNCTION + Real EntropyFromDensityTemperature(const Real rho, const Real temperature, + Real *lambda = nullptr) const { + return mpark::visit( + [&rho, &temperature, &lambda](const auto &eos) { + return eos.EntropyFromDensityTemperature(rho, temperature, lambda); + }, + eos_); + } + + PORTABLE_INLINE_FUNCTION + Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const { + return mpark::visit( + [&rho, &sie, &lambda](const auto &eos) { + return eos.EntropyFromDensityInternalEnergy(rho, sie, lambda); + }, + eos_); + } PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(const Real rho, const Real temperature, @@ -424,6 +443,93 @@ class Variant { eos_); } + template + inline void + EntropyFromDensityTemperature(ConstRealIndexer &&rhos, ConstRealIndexer &&temperatures, + RealIndexer &&entropies, const int num) const { + NullIndexer lambdas{}; // Returns null pointer for every index + return EntropyFromDensityTemperature(rhos, temperatures, entropies, num, lambdas); + } + + template + inline void EntropyFromDensityTemperature(ConstRealIndexer &&rhos, + ConstRealIndexer &&temperatures, + RealIndexer &&entropies, const int num, + LambdaIndexer &&lambdas) const { + return mpark::visit( + [&rhos, &temperatures, &entropies, &num, &lambdas](const auto &eos) { + return eos.EntropyFromDensityTemperature(rhos, temperatures, entropies, num, + lambdas); + }, + eos_); + } + + template + inline void EntropyFromDensityTemperature(ConstRealIndexer &&rhos, + ConstRealIndexer &&temperatures, + RealIndexer &&entropies, Real *scratch, + const int num) const { + NullIndexer lambdas{}; // Returns null pointer for every index + return EntropyFromDensityTemperature(rhos, temperatures, entropies, scratch, num, + lambdas); + } + + template + inline void + EntropyFromDensityTemperature(ConstRealIndexer &&rhos, ConstRealIndexer &&temperatures, + RealIndexer &&entropies, Real *scratch, const int num, + LambdaIndexer &&lambdas) const { + return mpark::visit( + [&rhos, &temperatures, &entropies, &scratch, &num, &lambdas](const auto &eos) { + return eos.EntropyFromDensityTemperature(rhos, temperatures, entropies, scratch, + num, lambdas); + }, + eos_); + } + + template + inline void + EntropyFromDensityInternalEnergy(ConstRealIndexer &&rhos, ConstRealIndexer &&sies, + RealIndexer &&entropies, const int num) const { + NullIndexer lambdas{}; // Returns null pointer for every index + return EntropyFromDensityInternalEnergy(rhos, sies, entropies, num, lambdas); + } + + template + inline void EntropyFromDensityInternalEnergy(ConstRealIndexer &&rhos, + ConstRealIndexer &&sies, + RealIndexer &&entropies, const int num, + LambdaIndexer &&lambdas) const { + return mpark::visit( + [&rhos, &sies, &entropies, &num, &lambdas](const auto &eos) { + return eos.EntropyFromDensityInternalEnergy(rhos, sies, entropies, num, + lambdas); + }, + eos_); + } + + template + inline void EntropyFromDensityInternalEnergy(ConstRealIndexer &&rhos, + ConstRealIndexer &&sies, + RealIndexer &&entropies, Real *scratch, + const int num) const { + NullIndexer lambdas{}; // Returns null pointer for every index + return EntropyFromDensityInternalEnergy(rhos, sies, entropies, scratch, num, lambdas); + } + + template + inline void + EntropyFromDensityInternalEnergy(ConstRealIndexer &&rhos, ConstRealIndexer &&sies, + RealIndexer &&entropies, Real *scratch, const int num, + LambdaIndexer &&lambdas) const { + return mpark::visit( + [&rhos, &sies, &entropies, &scratch, &num, &lambdas](const auto &eos) { + return eos.EntropyFromDensityInternalEnergy(rhos, sies, entropies, scratch, num, + lambdas); + }, + eos_); + } + template inline void SpecificHeatFromDensityTemperature(ConstRealIndexer &&rhos, ConstRealIndexer &&temperatures, diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index 25dea5b8a9..5534ffaa8a 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -55,6 +55,8 @@ class UnitSystem : public EosBase> { using EosBase>::InternalEnergyFromDensityTemperature; using EosBase>::PressureFromDensityTemperature; using EosBase>::PressureFromDensityInternalEnergy; + using EosBase>::EntropyFromDensityTemperature; + using EosBase>::EntropyFromDensityInternalEnergy; using EosBase>::SpecificHeatFromDensityTemperature; using EosBase>::SpecificHeatFromDensityInternalEnergy; using EosBase>::BulkModulusFromDensityTemperature; @@ -79,10 +81,11 @@ class UnitSystem : public EosBase> { const Real sie_unit, const Real temp_unit) : t_(std::forward(t)), rho_unit_(rho_unit), sie_unit_(sie_unit), temp_unit_(temp_unit), press_unit_(rho_unit * sie_unit), + entropy_unit_(sie_unit / temp_unit), inv_rho_unit_(1 / rho_unit) // inverses computed to avoid division at runtime , inv_sie_unit_(1 / sie_unit), inv_temp_unit_(1 / temp_unit), - inv_press_unit_(1 / press_unit_), + inv_press_unit_(1 / press_unit_), inv_entropy_unit_(1 / entropy_unit_), inv_dpde_unit_(sie_unit / press_unit_) // thermo derivatives computed consistently , inv_dvdt_unit_(rho_unit * temp_unit) // TODO(JMM): Is this convention weird? @@ -129,6 +132,13 @@ class UnitSystem : public EosBase> { return inv_press_unit_ * P; } PORTABLE_FUNCTION + Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const { + const Real S = + t_.EntropyFromDensityInternalEnergy(rho * rho_unit_, sie * sie_unit_, lambda); + return inv_entropy_unit_ * S; + } + PORTABLE_FUNCTION Real SpecificHeatFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda = nullptr) const { const Real cv = t_.SpecificHeatFromDensityInternalEnergy(rho * rho_unit_, @@ -156,6 +166,12 @@ class UnitSystem : public EosBase> { t_.PressureFromDensityTemperature(rho * rho_unit_, temp * temp_unit_, lambda); return inv_press_unit_ * P; } + Real EntropyFromDensityTemperature(const Real rho, const Real temp, + Real *lambda = nullptr) const { + const Real S = + t_.EntropyFromDensityTemperature(rho * rho_unit_, temp * temp_unit_, lambda); + return inv_entropy_unit_ * S; + } PORTABLE_FUNCTION Real SpecificHeatFromDensityTemperature(const Real rho, const Real temp, Real *lambda = nullptr) const { @@ -257,8 +273,8 @@ class UnitSystem : public EosBase> { T t_; // JMM: The compiler deletes GetOnDevice if I make these const. So... they're not - Real rho_unit_, sie_unit_, temp_unit_, press_unit_; - Real inv_rho_unit_, inv_sie_unit_, inv_temp_unit_, inv_press_unit_; + Real rho_unit_, sie_unit_, temp_unit_, press_unit_, entropy_unit_; + Real inv_rho_unit_, inv_sie_unit_, inv_temp_unit_, inv_press_unit_, inv_entropy_unit_; Real inv_dpde_unit_, inv_dvdt_unit_, inv_dpdr_unit_, inv_dtdr_unit_, inv_dtde_unit_; Real inv_cv_unit_, inv_bmod_unit_; }; diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 427a8aba49..6b90b08a1c 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -143,6 +143,11 @@ class BilinearRampEOS : public EosBase> { return p_eos < p_ramp ? p_ramp : p_eos; } PORTABLE_FUNCTION + Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const { + return t_.EntropyFromDensityInternalEnergy(rho, sie, lambda); + } + PORTABLE_FUNCTION Real SpecificHeatFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda = nullptr) const { return t_.SpecificHeatFromDensityInternalEnergy(rho, sie, lambda); @@ -169,6 +174,11 @@ class BilinearRampEOS : public EosBase> { return p_eos < p_ramp ? p_ramp : p_eos; } PORTABLE_FUNCTION + Real EntropyFromDensityTemperature(const Real rho, const Real temperature, + Real *lambda = nullptr) const { + return t_.EntropyFromDensityTemperature(rho, temperature, lambda); + } + PORTABLE_FUNCTION Real SpecificHeatFromDensityTemperature(const Real rho, const Real temperature, Real *lambda = nullptr) const { return t_.SpecificHeatFromDensityTemperature(rho, temperature, lambda); diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp index 8df9da8201..cbfec776dc 100644 --- a/singularity-eos/eos/modifiers/relativistic_eos.hpp +++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp @@ -47,6 +47,8 @@ class RelativisticEOS : public EosBase> { using EosBase>::InternalEnergyFromDensityTemperature; using EosBase>::PressureFromDensityTemperature; using EosBase>::PressureFromDensityInternalEnergy; + using EosBase>::EntropyFromDensityTemperature; + using EosBase>::EntropyFromDensityInternalEnergy; using EosBase>::SpecificHeatFromDensityTemperature; using EosBase>::SpecificHeatFromDensityInternalEnergy; using EosBase>::BulkModulusFromDensityTemperature; @@ -92,6 +94,11 @@ class RelativisticEOS : public EosBase> { return t_.PressureFromDensityInternalEnergy(rho, sie, lambda); } PORTABLE_FUNCTION + Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const { + return t_.EntropyFromDensityInternalEnergy(rho, sie, lambda); + } + PORTABLE_FUNCTION Real SpecificHeatFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda = nullptr) const { return t_.SpecificHeatFromDensityInternalEnergy(rho, sie, lambda); @@ -115,6 +122,11 @@ class RelativisticEOS : public EosBase> { return t_.PressureFromDensityTemperature(rho, temperature, lambda); } PORTABLE_FUNCTION + Real EntropyFromDensityTemperature(const Real rho, const Real temperature, + Real *lambda = nullptr) const { + return t_.EntropyFromDensityTemperature(rho, temperature, lambda); + } + PORTABLE_FUNCTION Real SpecificHeatFromDensityTemperature(const Real rho, const Real temperature, Real *lambda = nullptr) const { return t_.SpecificHeatFromDensityTemperature(rho, temperature, lambda); diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp index b04122fc9f..15590039ac 100644 --- a/singularity-eos/eos/modifiers/scaled_eos.hpp +++ b/singularity-eos/eos/modifiers/scaled_eos.hpp @@ -47,6 +47,8 @@ class ScaledEOS : public EosBase> { using EosBase>::InternalEnergyFromDensityTemperature; using EosBase>::PressureFromDensityTemperature; using EosBase>::PressureFromDensityInternalEnergy; + using EosBase>::EntropyFromDensityTemperature; + using EosBase>::EntropyFromDensityInternalEnergy; using EosBase>::SpecificHeatFromDensityTemperature; using EosBase>::SpecificHeatFromDensityInternalEnergy; using EosBase>::BulkModulusFromDensityTemperature; @@ -90,6 +92,11 @@ class ScaledEOS : public EosBase> { Real *lambda = nullptr) const { return t_.PressureFromDensityInternalEnergy(scale_ * rho, inv_scale_ * sie, lambda); } + Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const { + return scale_ * + t_.EntropyFromDensityInternalEnergy(scale_ * rho, inv_scale_ * sie, lambda); + } PORTABLE_FUNCTION Real SpecificHeatFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda = nullptr) const { @@ -114,6 +121,11 @@ class ScaledEOS : public EosBase> { return t_.PressureFromDensityTemperature(scale_ * rho, temperature, lambda); } PORTABLE_FUNCTION + Real EntropyFromDensityTemperature(const Real rho, const Real temperature, + Real *lambda = nullptr) const { + return scale_ * t_.EntropyFromDensityTemperature(scale_ * rho, temperature, lambda); + } + PORTABLE_FUNCTION Real SpecificHeatFromDensityTemperature(const Real rho, const Real temperature, Real *lambda = nullptr) const { return t_.SpecificHeatFromDensityTemperature(scale_ * rho, temperature, lambda); diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index 32666c2280..363db6c60e 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -47,6 +47,8 @@ class ShiftedEOS : public EosBase> { using EosBase>::InternalEnergyFromDensityTemperature; using EosBase>::PressureFromDensityTemperature; using EosBase>::PressureFromDensityInternalEnergy; + using EosBase>::EntropyFromDensityTemperature; + using EosBase>::EntropyFromDensityInternalEnergy; using EosBase>::SpecificHeatFromDensityTemperature; using EosBase>::SpecificHeatFromDensityInternalEnergy; using EosBase>::BulkModulusFromDensityTemperature; @@ -89,6 +91,11 @@ class ShiftedEOS : public EosBase> { return t_.PressureFromDensityInternalEnergy(rho, sie - shift_, lambda); } PORTABLE_FUNCTION + Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const { + return t_.EntropyFromDensityInternalEnergy(rho, sie - shift_, lambda); + } + PORTABLE_FUNCTION Real SpecificHeatFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda = nullptr) const { return t_.SpecificHeatFromDensityInternalEnergy(rho, sie - shift_, lambda); @@ -109,6 +116,11 @@ class ShiftedEOS : public EosBase> { return t_.PressureFromDensityTemperature(rho, temperature, lambda); } PORTABLE_FUNCTION + Real EntropyFromDensityTemperature(const Real rho, const Real temperature, + Real *lambda = nullptr) const { + return t_.EntropyFromDensityTemperature(rho, temperature, lambda); + } + PORTABLE_FUNCTION Real SpecificHeatFromDensityTemperature(const Real rho, const Real temperature, Real *lambda = nullptr) const { return t_.SpecificHeatFromDensityTemperature(rho, temperature, lambda); diff --git a/spack-repo/packages/singularity-eos/package.py b/spack-repo/packages/singularity-eos/package.py index 20f0df5323..d8045c1f96 100644 --- a/spack-repo/packages/singularity-eos/package.py +++ b/spack-repo/packages/singularity-eos/package.py @@ -64,7 +64,7 @@ class SingularityEos(CMakePackage, CudaPackage): depends_on("eospac", when="+eospac") depends_on("spiner") - depends_on("ports-of-call@1.4.1:") + depends_on("ports-of-call@1.4.2:") depends_on("spiner +kokkos", when="+kokkos") depends_on("mpark-variant") diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index 0a34072637..c91ea9bd37 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -32,6 +32,37 @@ PORTABLE_INLINE_FUNCTION Real QuadFormulaMinus(Real a, Real b, Real c) { return (-b - std::sqrt(b * b - 4 * a * c)) / (2 * a); } +// Only run exception checking when we aren't offloading to the GPUs +#ifdef PORTABILITY_STRATEGY_NONE +SCENARIO("Gruneisen EOS entropy is disabled", "[GruneisenEOS][Entropy]") { + GIVEN("Parameters for a Gruneisen EOS") { + // Unit conversions + constexpr Real cm = 1.; + constexpr Real us = 1e-06; + constexpr Real Mbcc_per_g = 1e12; + // Gruneisen parameters for copper + constexpr Real C0 = 0.394 * cm / us; + constexpr Real S1 = 1.489; + constexpr Real S2 = 0.; + constexpr Real S3 = 0.; + constexpr Real Gamma0 = 2.02; + constexpr Real b = 0.47; + constexpr Real rho0 = 8.93; + constexpr Real T0 = 298.; + constexpr Real P0 = 0.; + constexpr Real Cv = 0.383e-05 * Mbcc_per_g; + // Create the EOS + EOS host_eos = Gruneisen(C0, S1, S2, S3, Gamma0, b, rho0, T0, P0, Cv); + EOS eos = host_eos.GetOnDevice(); + THEN("A call to the entropy should throw an exception") { + using Catch::Matchers::Contains; + auto msg_matcher = Contains("Entropy is not enabled"); + REQUIRE_THROWS_WITH(eos.EntropyFromDensityTemperature(1.0, 1.0), msg_matcher); + } + } +} +#endif + SCENARIO("Gruneisen EOS", "[VectorEOS][GruneisenEOS]") { GIVEN("Parameters for a Gruneisen EOS") { // Unit conversions diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index 59c85ba0b1..a02c66230f 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -467,6 +467,45 @@ SCENARIO("EOS Unit System", "[EOSBuilder][UnitSystem][IdealGas]") { } } +SCENARIO("Ideal gas entropy", "[IdealGas][Entropy]") { + GIVEN("Parameters for an ideal gas with entropy reference states") { + // Create ideal gas EOS ojbect + constexpr Real Cv = 5.0; + constexpr Real gm1 = 0.4; + constexpr Real EntropyT0 = 100; + constexpr Real EntropyRho0 = 1e-03; + EOS host_eos = IdealGas(gm1, Cv, EntropyT0, EntropyRho0); + EOS eos = host_eos.GetOnDevice(); + THEN("The entropy at the reference state should be zero") { + auto entropy = eos.EntropyFromDensityTemperature(EntropyRho0, EntropyT0); + INFO("Entropy should be zero but it is " << entropy); + CHECK(isClose(entropy, 0.0, 1.e-14)); + } + GIVEN("A state at the reference temperature and a density whose cube root is the " + "reference density") { + constexpr Real T = EntropyT0; + constexpr Real rho = 0.1; // rho**3 = EntropyRho0 + THEN("The entropy should be 2. / 3. * gm1 * Cv * log(EntropyRho0)") { + const Real entropy_true = 2. / 3. * gm1 * Cv * log(EntropyRho0); + auto entropy = eos.EntropyFromDensityTemperature(rho, T); + INFO("Entropy: " << entropy << " True entropy: " << entropy_true); + CHECK(isClose(entropy, entropy_true, 1e-12)); + } + } + GIVEN("A state at the reference density and a temperature whose square is the " + "reference temperature") { + constexpr Real T = 10; // T**2 = EntropyT0 + constexpr Real rho = EntropyRho0; + THEN("The entropy should be -1. / 2. * Cv * log(EntropyT0)") { + const Real entropy_true = -1. / 2. * Cv * log(EntropyT0); + auto entropy = eos.EntropyFromDensityTemperature(rho, T); + INFO("Entropy: " << entropy << " True entropy: " << entropy_true); + CHECK(isClose(entropy, entropy_true, 1e-12)); + } + } + } +} + SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") { GIVEN("Parameters for an ideal gas") { @@ -517,15 +556,38 @@ SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") { constexpr std::array heatcapacity_true{Cv, Cv, Cv}; constexpr std::array gruneisen_true{gm1, gm1, gm1}; + // Gold standard entropy doesn't produce round numbers so we need to + // calculate it from the device views so this requires a bit more work +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View v_entropy_true("True Entropy"); +#else + std::array entropy_true; + Real *v_entropy_true = entropy_true.data(); +#endif + constexpr Real P0 = 1e6; // microbar + constexpr Real T0 = 293; // K + constexpr Real rho0 = P0 / (gm1 * Cv * T0); // g/cm^3 + portableFor( + "Calculate true entropy", 0, num, PORTABLE_LAMBDA(const int i) { + v_entropy_true[i] = + Cv * log(v_energy[i] / Cv / T0) + gm1 * Cv * log(rho0 / v_density[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + auto entropy_true = Kokkos::create_mirror_view(v_entropy_true); + Kokkos::deep_copy(entropy_true, v_entropy_true); +#endif + #ifdef PORTABILITY_STRATEGY_KOKKOS // Create device views for outputs and mirror those views on the host Kokkos::View v_temperature("Temperature"); Kokkos::View v_pressure("Pressure"); + Kokkos::View v_entropy("Entropy"); Kokkos::View v_heatcapacity("Cv"); Kokkos::View v_bulkmodulus("bmod"); Kokkos::View v_gruneisen("Gamma"); auto h_temperature = Kokkos::create_mirror_view(v_temperature); auto h_pressure = Kokkos::create_mirror_view(v_pressure); + auto h_entropy = Kokkos::create_mirror_view(v_entropy); auto h_heatcapacity = Kokkos::create_mirror_view(v_heatcapacity); auto h_bulkmodulus = Kokkos::create_mirror_view(v_bulkmodulus); auto h_gruneisen = Kokkos::create_mirror_view(v_gruneisen); @@ -534,12 +596,14 @@ SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") { // will be passed to the functions in place of the Kokkos views std::array h_temperature; std::array h_pressure; + std::array h_entropy; std::array h_heatcapacity; std::array h_bulkmodulus; std::array h_gruneisen; // Just alias the existing pointers auto v_temperature = h_temperature.data(); auto v_pressure = h_pressure.data(); + auto v_entropy = h_entropy.data(); auto v_heatcapacity = h_heatcapacity.data(); auto v_bulkmodulus = h_bulkmodulus.data(); auto v_gruneisen = h_gruneisen.data(); @@ -570,6 +634,18 @@ SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") { } } + WHEN("An S(rho, e) lookup is performed") { + eos.EntropyFromDensityInternalEnergy(v_density, v_energy, v_entropy, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_entropy, v_entropy); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned S(rho, e) should be equal to the true entropy") { + array_compare(num, density, energy, h_entropy, entropy_true, "Density", + "Energy"); + } + } + WHEN("A C_v(rho, e) lookup is performed") { eos.SpecificHeatFromDensityInternalEnergy(v_density, v_energy, v_heatcapacity, num); @@ -650,15 +726,38 @@ SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") { constexpr std::array heatcapacity_true{Cv, Cv, Cv}; constexpr std::array gruneisen_true{gm1, gm1, gm1}; + // Gold standard entropy doesn't produce round numbers so we need to + // calculate it from the device views so this requires a bit more work +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View v_entropy_true("True Entropy"); +#else + std::array entropy_true; + Real *v_entropy_true = entropy_true.data(); +#endif + constexpr Real P0 = 1e6; // microbar + constexpr Real T0 = 293; // K + constexpr Real rho0 = P0 / (gm1 * Cv * T0); // g/cm^3 + portableFor( + "Calculate true entropy", 0, num, PORTABLE_LAMBDA(const int i) { + v_entropy_true[i] = + Cv * log(v_temperature[i] / T0) + gm1 * Cv * log(rho0 / v_density[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + auto entropy_true = Kokkos::create_mirror_view(v_entropy_true); + Kokkos::deep_copy(entropy_true, v_entropy_true); +#endif + #ifdef PORTABILITY_STRATEGY_KOKKOS // Create device views for outputs and mirror those views on the host Kokkos::View v_energy("Energy"); Kokkos::View v_pressure("Pressure"); + Kokkos::View v_entropy("Entropy"); Kokkos::View v_heatcapacity("Cv"); Kokkos::View v_bulkmodulus("bmod"); Kokkos::View v_gruneisen("Gamma"); auto h_energy = Kokkos::create_mirror_view(v_energy); auto h_pressure = Kokkos::create_mirror_view(v_pressure); + auto h_entropy = Kokkos::create_mirror_view(v_entropy); auto h_heatcapacity = Kokkos::create_mirror_view(v_heatcapacity); auto h_bulkmodulus = Kokkos::create_mirror_view(v_bulkmodulus); auto h_gruneisen = Kokkos::create_mirror_view(v_gruneisen); @@ -667,12 +766,14 @@ SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") { // will be passed to the functions in place of the Kokkos views std::array h_energy; std::array h_pressure; + std::array h_entropy; std::array h_heatcapacity; std::array h_bulkmodulus; std::array h_gruneisen; // Just alias the existing pointers auto v_energy = h_energy.data(); auto v_pressure = h_pressure.data(); + auto v_entropy = h_entropy.data(); auto v_heatcapacity = h_heatcapacity.data(); auto v_bulkmodulus = h_bulkmodulus.data(); auto v_gruneisen = h_gruneisen.data(); @@ -702,6 +803,18 @@ SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") { } } + WHEN("An S(rho, T) lookup is performed") { + eos.EntropyFromDensityTemperature(v_density, v_temperature, v_entropy, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_entropy, v_entropy); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned S(rho, T) should be equal to the true entropy") { + array_compare(num, density, temperature, h_entropy, entropy_true, "Density", + "Temperature"); + } + } + WHEN("A C_v(rho, T) lookup is performed") { eos.SpecificHeatFromDensityTemperature(v_density, v_temperature, v_heatcapacity, num);