From 45ee5d634197fe3c30108bce091e9048a2ee14b5 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 09:34:26 -0700 Subject: [PATCH 01/44] Add sublime text and Mac OSX filesystem things to the gitignore --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) 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* + From 2228efbc67c3be51f0d5f9472292d383c8b70848 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 09:35:49 -0700 Subject: [PATCH 02/44] Remove commented code, add vector entropy functions, and add a default function for entropy not being enabled --- singularity-eos/eos/eos_base.hpp | 103 ++++++++++++------------------- 1 file changed, 39 insertions(+), 64 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index d94ef5ca0f..6e4b4b8726 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -18,6 +18,7 @@ #include #include +#include namespace singularity { namespace mfuncname { @@ -50,7 +51,9 @@ namespace eos_base { using EosBase::MinimumDensity; \ using EosBase::MinimumTemperature; \ using EosBase::PTofRE; \ - using EosBase::FillEos; + using EosBase::FillEos; \ + using EosBase::EntropyFromDensityTemperature; \ + using EosBase::EntropyFromDensityEnergy; /* This is a CRTP that allows for static inheritance so that default behavior for @@ -121,6 +124,34 @@ class EosBase { }); } 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 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 SpecificHeatFromDensityTemperature(ConstRealIndexer &&rhos, ConstRealIndexer &&temperatures, RealIndexer &&cvs, const int num, @@ -264,69 +295,13 @@ 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 + [[noreturn]] PORTABLE_FORCE_INLINE_FUNCTION + void EntropyIsNotEnabled() const { + PORTABLE_ALWAYS_THROW_OR_ABORT(std::string("Entropy is not enabled for the '") + + std::string(typeid(CRTP).name()) + + std::string("' EOS")); + } }; } // namespace eos_base } // namespace singularity From f408ec51206aa31bdc310e152cf5549996b7ab98 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 09:36:21 -0700 Subject: [PATCH 03/44] Add entropy functions to variant --- singularity-eos/eos/eos_variant.hpp | 61 +++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index 0fbe1b8312..cdbbeec16f 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, @@ -332,6 +351,48 @@ 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 + 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 SpecificHeatFromDensityTemperature(ConstRealIndexer &&rhos, ConstRealIndexer &&temperatures, From a4441676309e6647586d84475647a1dc8ba0bd86 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 09:37:24 -0700 Subject: [PATCH 04/44] Add section on complete EOS and the beginnings of documentation for entropy --- doc/sphinx/src/models.rst | 176 +++++++++++++++++++++++++++++++++----- 1 file changed, 156 insertions(+), 20 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index ef5269b2d6..b65ee42796 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -7,7 +7,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,8 +22,8 @@ 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 @@ -59,8 +59,100 @@ 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 ---------------------- +Completing an Incomplete Equation of State and Consequences for 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. + +In general, 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, :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, the Mie Gruneisen +form tends to dominate most analytic EOS and it becomes imperitive to extend the +incomplete form to a complete equation of state. + +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 :ref:`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 :ref:`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 +:ref:`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. + +.. _Davis reactants: https://doi.org/10.1016/S0010-2180(99)00112-1 + +.. _Maxwell relations: https://en.wikipedia.org/wiki/Maxwell_relations + +.. _Menikoff Complete EOS: + https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-UR-16-21706 + +Again from :ref:`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 @@ -72,10 +164,14 @@ in the :doc:`EOS API `. | :math:`T(\rho, e)` | Temperature | Density, Internal Energy | +--------------------------+----------------------+ | | :math:`P(\rho, e)` | Pressure | | ++--------------------------+----------------------+ | +| :math:`S(\rho, e)` | Entropy | | +--------------------------+----------------------+--------------------------+ | :math:`e(\rho, T)` | Internal Energy | Density, Temperature | +--------------------------+----------------------+ | | :math:`P(\rho, T)` | Pressure | | ++--------------------------+----------------------+ | +| :math:`S(\rho, T)` | Entropy | | +--------------------------+----------------------+--------------------------+ | :math:`\rho(P, T)` | Density | Pressure, Temperature | +--------------------------+----------------------+ | @@ -99,11 +195,27 @@ 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). -Disambiguation -```````````````` - -Gruneisen Parameter -''''''''''''''''''' +Entropy availability +```````````````````` +Although the functionality to calculate an entropy is exposed for all equations +of state, the fact that many equations of state are based off of incomplete +equations of state means that the entropy for many EOS is not thermodynamically +consistent. While this is also true of the temperature, we have found that +thermodynamic consistency is often more important for calculations involving +entropy. + +As a result, we *only* expose a numerical value for the entropy when the +equation of state form is fully complete and thermodynamically constistent. We +should note the standard sesame *form* is generally thermodynamically +consistent even if specific data tables may not be. We do not *presently* +provide tools to determine the thermodynamic consistency, either globally or +locally, for a sesame EOS, but we would like to in the future. + +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 +259,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 +270,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,10 +293,25 @@ the form e = C_V T, -where quantities are defined in the :ref:`nomenclature ` section. +where quantities are defined in the :ref:`nomenclature ` section. -The settable parameters are the Gruneisen parameter and specific heat capacity. +Generally, the change in entropy for an ideal gas is given by +.. math:: + + \Delta S = C_V \ln\left(\frac{T}{T_0}\right) + R \ln\left(\frac{\rho_0} + {\rho}\right), + +where :math:`R` is the *specific* gas constant in units of erg/g-K. However, it +should be noted that this equation may be invalid for temperatures below the +reference temperature as the first term will tend towards negative infinity. +This is limitation of assuming a constant heat capacity since the third law of +thermodynamics states that the heat capacity should approach zero as the +temperature approaches zero, whereupon the entropy should be constant. + +The settable parameters for this EOS are the Gruneisen parameter and specific +heat capacity. 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 @@ -222,8 +351,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 +367,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 @@ -393,21 +526,24 @@ 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 +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. + +Using the fact that the heat capacity is defined by The :math:`e(\rho, P)` lookup is quite awkward, so the energy is more-conveniently cast in terms of termperature such that @@ -429,7 +565,7 @@ The Gruneisen parameter takes on a linear form such that Zy & \rho >= \rho_0 \end{cases} -where :math:`Z` and :math:`y` are dimensionless parameters. +where :math:`Z` is a dimensionless parameter and :math:`y = 1 - \rho0/\rho`. Finally, the pressure, energy, and temperature along the isentrope are given by From 2d615d2ed0db20914c0560703f30af4d79ce240c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 15:36:36 -0700 Subject: [PATCH 05/44] Finalize entropy documentation --- doc/sphinx/src/models.rst | 312 +++++++++++++++++++++++++------------- 1 file changed, 207 insertions(+), 105 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index adfa91eff1..14895c9337 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 =========== @@ -26,11 +42,9 @@ 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 ```````````````````````` @@ -76,10 +90,10 @@ variables*, i.e. .. math:: - e = e(\rho, S) \\ - h = h(P, S) \\ - F = F(\rho, T) \\ - G = G(P, T) \\ + e =& e(\rho, S) \\ + h =& h(P, S) \\ + F =& F(\rho, T) \\ + G =& G(P, T) \\ where all potentials are specific, :math:`e` is again the internal energy, :math:`h` is the enthalpy, :math:`F` is the Helmholtz free energy, and :math:`G` @@ -102,11 +116,11 @@ specific requirements placed on the functional form of the heat capacity described below. -The :ref:`Maxwell relations ` are the consequence of the +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 :ref:`here -`) to produce the relationship, +can be extended to third derivatives (as Davis does in his Appendix B `here +`_) to produce the relationship, .. math:: @@ -114,7 +128,7 @@ can be extended to third derivatives (as Davis does in his Appendix B :ref:`here \left(\frac{\partial \Gamma}{\partial S}\right)_V. This is often referred to as a "compatibility condition" (see also -:ref:`Menikoff `) and provides an important connection +`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 @@ -122,14 +136,7 @@ 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. -.. _Davis reactants: https://doi.org/10.1016/S0010-2180(99)00112-1 - -.. _Maxwell relations: https://en.wikipedia.org/wiki/Maxwell_relations - -.. _Menikoff Complete EOS: - https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-UR-16-21706 - -Again from :ref:`Menikoff `, it can be shown that even +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 @@ -143,12 +150,12 @@ energy along the isentrope that passes through the reference temperature, and .. math:: - \phi(\rho) = \exp\left(\int\limits_{\rho_0}^{\rho} \rho' \Gamma - (\rho') \mathrm{d}\rho' \right). + \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 +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 @@ -161,55 +168,95 @@ in the :doc:`EOS API `. +--------------------------+----------------------+--------------------------+ | Function | Dependent Variable | Independent Variables | +==========================+======================+==========================+ -| :math:`T(\rho, e)` | Temperature | Density, Internal Energy | -+--------------------------+----------------------+ | +| :math:`T(\rho, e)` | Temperature | Density, Specific | ++--------------------------+----------------------+ Internal Energy | | :math:`P(\rho, e)` | Pressure | | +--------------------------+----------------------+ | -| :math:`S(\rho, e)` | Entropy | | +| :math:`S(\rho, e)` | Specific Entropy | | +--------------------------+----------------------+--------------------------+ -| :math:`e(\rho, T)` | Internal Energy | Density, Temperature | +| :math:`e(\rho, T)` | Specific Internal | Density, Temperature | +| | Energy | | +--------------------------+----------------------+ | | :math:`P(\rho, T)` | Pressure | | +--------------------------+----------------------+ | -| :math:`S(\rho, T)` | Entropy | | +| :math:`S(\rho, T)` | Specific Entropy | | +--------------------------+----------------------+--------------------------+ | :math:`\rho(P, T)` | Density | Pressure, Temperature | +--------------------------+----------------------+ | -| :math:`e(P, T)` | Internal Energy | | +| :math:`e(P, T)` | Specific Internal | | +| | Energy | | +--------------------------+----------------------+--------------------------+ | :math:`C_V(\rho, T)` | Constant Volume | Density, Temperature | +--------------------------+ Specific Heat +--------------------------+ -| :math:`C_V(\rho, e)` | Capacity | Density, Internal Energy | +| :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, Internal Energy | +| :math:`B_S(\rho, e)` | | Density, Specific | +| | | Internal Energy | +--------------------------+----------------------+--------------------------+ | :math:`\Gamma(\rho, T)` | Gruneisen Parameter | Density, Temperature | +--------------------------+ +--------------------------+ -| :math:`\Gamma(\rho, e)` | | Density, Internal Energy | +| :math:`\Gamma(\rho, e)` | | Density, Specific | +| | | Internal Energy | +--------------------------+----------------------+--------------------------+ 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). Entropy availability ```````````````````` -Although the functionality to calculate an entropy is exposed for all equations -of state, the fact that many equations of state are based off of incomplete -equations of state means that the entropy for many EOS is not thermodynamically -consistent. While this is also true of the temperature, we have found that -thermodynamic consistency is often more important for calculations involving -entropy. - -As a result, we *only* expose a numerical value for the entropy when the -equation of state form is fully complete and thermodynamically constistent. We -should note the standard sesame *form* is generally thermodynamically -consistent even if specific data tables may not be. We do not *presently* -provide tools to determine the thermodynamic consistency, either globally or -locally, for a sesame EOS, but we would like to in the future. +For an arbitrary equation of state, a change in entropy in terms of temperature +and volume is given by + +.. 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 ```````````````````````````` @@ -282,6 +329,12 @@ Implemented EOS models Ideal Gas ````````` +.. note:: + + While entropy is definied for this EOS, it diverges to negative infinity at + absolute zero. Care should be taken using values significantly below that of + the reference state. + The ideal gas (aka perfect or gamma-law gas) model in ``singularity-eos`` takes the form @@ -294,39 +347,50 @@ the form e = C_V T, where quantities are defined in the :ref:`nomenclature ` section. +Information and 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 + +.. math:: -Generally, the change in entropy for an ideal gas is given by + \Gamma = \gamma - 1 + +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. + +The entropy for an ideal gas is given by .. math:: - \Delta S = C_V \ln\left(\frac{T}{T_0}\right) + R \ln\left(\frac{\rho_0} + S = C_V \ln\left(\frac{T}{T_0}\right) + \Gamma C_V \ln\left(\frac{\rho_0} {\rho}\right), -where :math:`R` is the *specific* gas constant in units of erg/g-K. However, it -should be noted that this equation may be invalid for temperatures below the -reference temperature as the first term will tend towards negative infinity. -This is limitation of assuming a constant heat capacity since the third law of -thermodynamics states that the heat capacity should approach zero as the -temperature approaches zero, whereupon the entropy should be constant. +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 -The settable parameters for this EOS are the Gruneisen parameter and specific -heat capacity. +.. math:: -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. + \rho_0 = \frac{P_0}{\Gamma C_V T_0}, -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 +where :math:`P_0 = 1` bar. -.. math:: +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. - \Gamma = \gamma - 1 +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 @@ -334,11 +398,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 @@ -418,7 +498,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 @@ -459,17 +540,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:: @@ -483,19 +564,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:: @@ -505,10 +587,12 @@ 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:: @@ -526,13 +610,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 @@ -579,12 +668,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 @@ -604,9 +694,12 @@ Information and Nomenclature>` section. Davis Reactants 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 +.. 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:: @@ -692,6 +785,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 @@ -767,6 +863,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: @@ -978,6 +1077,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 From d2110d6d1e62be5d95580123e7b5c3d6b4a3adb6 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 15:37:14 -0700 Subject: [PATCH 06/44] Add entropy to the ideal gas --- singularity-eos/eos/eos_ideal.hpp | 41 +++++++++++++++++++++++++++++-- 1 file changed, 39 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp index 0c6ed01dd6..c2a6221a62 100644 --- a/singularity-eos/eos/eos_ideal.hpp +++ b/singularity-eos/eos/eos_ideal.hpp @@ -19,6 +19,7 @@ #include #include #include +#include // Ports-of-call #include @@ -26,6 +27,7 @@ // Base stuff #include #include +#include #include #define MYMAX(a, b) a > b ? a : b @@ -40,13 +42,36 @@ 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) + { + checkInputs(); + } - IdealGas GetOnDevice() { return *this; } + 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); + PORTABLE_ALWAYS_REQUIRE(_EntropyT0 <= 0); + PORTABLE_ALWAYS_REQUIRE(_EntropyRho0 <= 0); + } PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const { return MYMAX(0.0, _Cv * temperature); @@ -59,6 +84,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, density)); + } + 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, sie, lambda); + } PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const { return _Cv; @@ -128,6 +163,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 From 5b4be7879ea39e02d0feee03ff4c970bde4b8fae Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 16:10:54 -0700 Subject: [PATCH 07/44] Add EntropyIsNotEnabled to default using statements --- singularity-eos/eos/eos_base.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 6e4b4b8726..11ed54f8ad 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -53,7 +53,8 @@ namespace eos_base { using EosBase::PTofRE; \ using EosBase::FillEos; \ using EosBase::EntropyFromDensityTemperature; \ - using EosBase::EntropyFromDensityEnergy; + using EosBase::EntropyFromDensityEnergy; \ + using EosBase::EntropyIsNotEnabled; /* This is a CRTP that allows for static inheritance so that default behavior for From a561bd064e9aaf1d9c700f99f7b3ab8b47623c5e Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 16:11:30 -0700 Subject: [PATCH 08/44] Disable entropy in other EOS for the time being --- singularity-eos/eos/eos_davis.hpp | 10 +++++++ singularity-eos/eos/eos_eospac.hpp | 15 +++++++++++ singularity-eos/eos/eos_gruneisen.hpp | 16 +++++++++++ singularity-eos/eos/eos_spiner.hpp | 38 +++++++++++++++++++++++++++ 4 files changed, 79 insertions(+) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index b3291d2bc2..44c4c91b41 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(); + return 1.0 + } + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + EntropyIsNotEnabled(); + return 1.0 + } PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( const Real rho, const Real temp, Real *lambda = nullptr) const { return SpecificHeatFromDensityInternalEnergy( diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp index c629a033ac..f55f8c75ce 100644 --- a/singularity-eos/eos/eos_eospac.hpp +++ b/singularity-eos/eos/eos_eospac.hpp @@ -81,6 +81,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; @@ -89,6 +91,7 @@ class EOSPAC : public EosBase { using EosBase::GruneisenParamFromDensityInternalEnergy; using EosBase::PTofRE; using EosBase::FillEos; + using EosBase::EntropyIsNotEnabled; // TODO (JHP): Change EOSPAC vector implementations to be more performant // template @@ -273,6 +276,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 = nullptr) const { + EntropyIsNotEnabled(); + 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, @@ -386,6 +395,12 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::PressureFromDensityInternalEnergy( Real temp = TemperatureFromDensityInternalEnergy(rho, sie, lambda); return PressureFromDensityTemperature(rho, temp, lambda); } +PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) 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 f6de2f16c9..7c521651b8 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; @@ -306,6 +310,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(); + return 1.0; +} PORTABLE_INLINE_FUNCTION Real Gruneisen::BulkModulusFromDensityInternalEnergy( const Real rho_in, const Real sie, Real *lambda) const { using namespace gruneisen_utils; @@ -328,6 +338,12 @@ PORTABLE_INLINE_FUNCTION Real Gruneisen::PressureFromDensityTemperature( return PressureFromDensityInternalEnergy( rho, InternalEnergyFromDensityTemperature(rho, temp)); } +PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( + const Real rho_in, const Real sie, Real *lambda = nullptr) 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_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index 6ffd185d4c..1345294612 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -85,6 +85,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 +110,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 @@ -304,6 +311,7 @@ class SpinerEOSDependsRhoSie : public EosBase { using EosBase::GruneisenParamFromDensityInternalEnergy; using EosBase::PTofRE; using EosBase::FillEos; + using EosBase::EntropyIsNotEnabled; PORTABLE_INLINE_FUNCTION SpinerEOSDependsRhoSie() : memoryStatus_(DataStatus::Deallocated) {} @@ -924,6 +932,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(); + return 1.0; +} +PORTABLE_INLINE_FUNCTION +Real SpinerEOSDependsRhoT::EntropyFromDensityInternalEnergy(const Real rho, + const Real sie, + Real *lambda) const { + EntropyIsNotEnabled(); + return 1.0; +} + PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::SpecificHeatFromDensityTemperature(const Real rho, const Real temperature, @@ -1640,6 +1663,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(); + return 1.0; +} +PORTABLE_INLINE_FUNCTION +Real SpinerEOSDependsRhoSie::EntropyFromDensityInternalEnergy(const Real rho, + const Real sie, + Real *lambda) const { + EntropyIsNotEnabled(); + return 1.0; +} + PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoSie::SpecificHeatFromDensityTemperature(const Real rho, const Real T, From 7c98861017a3faf37cd0fe5c58b40d956e361567 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 17:25:41 -0700 Subject: [PATCH 09/44] Add entropy to unitsystem modifier --- .../eos/modifiers/eos_unitsystem.hpp | 21 +++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index 74e8174efb..5cc9d9f643 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,12 @@ 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_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 +133,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 +167,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 { @@ -249,8 +266,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_; }; From be3e6585101da6d44abf2d552a26a244bb4fc568 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 19:09:05 -0700 Subject: [PATCH 10/44] Add Helmholtz explanation for scaling ratios that leads to entropy scaling --- doc/sphinx/src/modifiers.rst | 39 +++++++++++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) 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 From 44a9529d04036bdd40c198fcf449265c239c8524 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 19:10:12 -0700 Subject: [PATCH 11/44] Add basic, unmodified, entropy lookups for ramps and relativistic modifiers --- singularity-eos/eos/modifiers/ramps_eos.hpp | 10 ++++++++++ singularity-eos/eos/modifiers/relativistic_eos.hpp | 10 ++++++++++ 2 files changed, 20 insertions(+) diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 567dc099ed..5b8bef549a 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 311d5af42d..83e07728df 100644 --- a/singularity-eos/eos/modifiers/relativistic_eos.hpp +++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp @@ -92,6 +92,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 +120,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); From 234a34c46ec62bee2f686655f01656c62db2b746 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 19:10:39 -0700 Subject: [PATCH 12/44] Add entropy lookups with shifted energies --- singularity-eos/eos/modifiers/shifted_eos.hpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index 1e27453bf2..6b51906159 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -89,6 +89,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 +114,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); From 3d94fc9ead8d00ca34ca7e9f2d9797fa0b8f2626 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 19:10:59 -0700 Subject: [PATCH 13/44] Apply scaling ratio to entropy output --- singularity-eos/eos/modifiers/scaled_eos.hpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp index 74427cd9ab..d054999a75 100644 --- a/singularity-eos/eos/modifiers/scaled_eos.hpp +++ b/singularity-eos/eos/modifiers/scaled_eos.hpp @@ -90,6 +90,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 +119,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); From a5cdf7cb335788834b78c3355dabfbdd36c9fdd2 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 19:25:32 -0700 Subject: [PATCH 14/44] Move note and add one to Vinet on divergent entropy --- doc/sphinx/src/models.rst | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 14895c9337..65a4557a5a 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -329,12 +329,6 @@ Implemented EOS models Ideal Gas ````````` -.. note:: - - While entropy is definied for this EOS, it diverges to negative infinity at - absolute zero. Care should be taken using values significantly below that of - the reference state. - The ideal gas (aka perfect or gamma-law gas) model in ``singularity-eos`` takes the form @@ -371,6 +365,12 @@ The entropy for an ideal gas is given by 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 @@ -599,8 +599,13 @@ coefficients :math:`d_n`, :math:`n\geq 2`, by 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 From eff8b4df00e5df7cecf4a3a999517573b094de56 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 19:34:13 -0700 Subject: [PATCH 15/44] Reorganize changelog --- CHANGELOG.md | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a44d158a49..8adcdd052b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,24 +2,18 @@ ## Current develop -### Fixed (Repair, bugs, ect) -- [[PR232]](https://github.com/lanl/singularity-eos/pull/228) Fixed uninitialized cmake path variables - -### Fixed (issue #227) -- [[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 - -### Added (entropy calculation for stellar collapse eos) -- [[PR226]](https://github.com/lanl/singularity-eos/pull/226) added entropy interpolation to stellar collapse eos - ### Fixed (Repair bugs, etc) - [[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/...) - [[PR202]](https://github.com/lanl/singularity-eos/pull/202) added the Vinet analytical EOS wth test cases and documentation. - [[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 +- [[PR226]](https://github.com/lanl/singularity-eos/pull/226) added entropy interpolation to stellar collapse 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 ### Changed (changing behavior/API/variables/...) - [[PR223]](https://github.com/lanl/singularity-eos/pull/223) Update ports-of-call and add portable error handling From 25c8d46016c6d29998d96b75b2d028851d6bb3fc Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 19:36:54 -0700 Subject: [PATCH 16/44] Update changelog for entropy additions --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8adcdd052b..e05a6ab6f5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,10 +14,12 @@ - [[PR214]](https://github.com/lanl/singularity-eos/pull/214) added documentation about adding a new EOS - [[PR226]](https://github.com/lanl/singularity-eos/pull/226) added entropy interpolation to stellar collapse 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 for the ideal gas and modifiers and an error for EOS where entropy is not implemented yet ### Changed (changing behavior/API/variables/...) - [[PR223]](https://github.com/lanl/singularity-eos/pull/223) Update ports-of-call and add portable 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 From 3b2ffdccd78c008b7aacc146e220c69c04ac3616 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Feb 2023 19:40:39 -0700 Subject: [PATCH 17/44] Whoops... cmath --- singularity-eos/eos/eos_ideal.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp index c2a6221a62..eb7ff976ef 100644 --- a/singularity-eos/eos/eos_ideal.hpp +++ b/singularity-eos/eos/eos_ideal.hpp @@ -19,7 +19,7 @@ #include #include #include -#include +#include // Ports-of-call #include From 05e7bf992263763447f5705709aac398f3b46706 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 17 Feb 2023 17:49:55 -0700 Subject: [PATCH 18/44] Fix typos and compile-time bugs --- singularity-eos/eos/eos_base.hpp | 4 ++-- singularity-eos/eos/eos_davis.hpp | 4 ++-- singularity-eos/eos/eos_eospac.hpp | 12 ++++++++---- singularity-eos/eos/eos_gruneisen.hpp | 4 ++-- singularity-eos/eos/eos_ideal.hpp | 10 +++++----- singularity-eos/eos/eos_spiner.hpp | 6 ++++++ singularity-eos/eos/modifiers/eos_unitsystem.hpp | 8 ++++---- 7 files changed, 29 insertions(+), 19 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 11ed54f8ad..28c71b0660 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -53,7 +53,7 @@ namespace eos_base { using EosBase::PTofRE; \ using EosBase::FillEos; \ using EosBase::EntropyFromDensityTemperature; \ - using EosBase::EntropyFromDensityEnergy; \ + using EosBase::EntropyFromDensityInternalEnergy; \ using EosBase::EntropyIsNotEnabled; /* @@ -297,7 +297,7 @@ class EosBase { Real RhoPmin(const Real temp) const { return 0.0; } // Default entropy behavior is to return an error - [[noreturn]] PORTABLE_FORCE_INLINE_FUNCTION + [[noreturn]] PORTABLE_FORCEINLINE_FUNCTION void EntropyIsNotEnabled() const { PORTABLE_ALWAYS_THROW_OR_ABORT(std::string("Entropy is not enabled for the '") + std::string(typeid(CRTP).name()) + diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 44c4c91b41..4382cb91cd 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -63,12 +63,12 @@ class DavisReactants : public EosBase { PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const { EntropyIsNotEnabled(); - return 1.0 + return 1.0; } PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda = nullptr) const { EntropyIsNotEnabled(); - return 1.0 + return 1.0; } PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( const Real rho, const Real temp, Real *lambda = nullptr) const { diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp index f55f8c75ce..a45ee0f829 100644 --- a/singularity-eos/eos/eos_eospac.hpp +++ b/singularity-eos/eos/eos_eospac.hpp @@ -50,6 +50,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( @@ -277,9 +281,9 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::PressureFromDensityTemperature(const Real } PORTABLE_INLINE_FUNCTION Real EOSPAC::EntropyFromDensityTemperature( - const Real rho, const Real temperature, Real *lambda = nullptr) const { + const Real rho, const Real temperature, Real *lambda) const { EntropyIsNotEnabled(); - return 1.0 + return 1.0; } PORTABLE_INLINE_FUNCTION void EOSPAC::FillEos(Real &rho, Real &temp, Real &sie, @@ -395,8 +399,8 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::PressureFromDensityInternalEnergy( Real temp = TemperatureFromDensityInternalEnergy(rho, sie, lambda); return PressureFromDensityTemperature(rho, temp, lambda); } -PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( - const Real rho, const Real sie, Real *lambda = nullptr) const { +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); diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index 7c521651b8..800c004cdc 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -338,8 +338,8 @@ PORTABLE_INLINE_FUNCTION Real Gruneisen::PressureFromDensityTemperature( return PressureFromDensityInternalEnergy( rho, InternalEnergyFromDensityTemperature(rho, temp)); } -PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( - const Real rho_in, const Real sie, Real *lambda = nullptr) const { +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); diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp index eb7ff976ef..1b196f3a60 100644 --- a/singularity-eos/eos/eos_ideal.hpp +++ b/singularity-eos/eos/eos_ideal.hpp @@ -54,7 +54,7 @@ class IdealGas : public EosBase { _dvdt0(1. / (_rho0 * _T0)), _EntropyT0(EntropyT0), _EntropyRho0(EntropyRho0) { - checkInputs(); + checkParams(); } IdealGas GetOnDevice() { @@ -68,9 +68,9 @@ class IdealGas : public EosBase { // 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); - PORTABLE_ALWAYS_REQUIRE(_EntropyT0 <= 0); - PORTABLE_ALWAYS_REQUIRE(_EntropyRho0 <= 0); + 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 { @@ -87,7 +87,7 @@ class IdealGas : public EosBase { 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, density)); + _gm1 * _Cv * log(robust::ratio(_EntropyRho0, rho)); } PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda = nullptr) const { diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index 1345294612..16fab3a15b 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -337,6 +337,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 diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index 5cc9d9f643..590a8e8ea8 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -81,12 +81,12 @@ 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), + 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_entropy_unit(1 / entropy_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? @@ -266,8 +266,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_, entropy_unit; - Real inv_rho_unit_, inv_sie_unit_, inv_temp_unit_, inv_press_unit_, inv_entropy_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_; }; From 6d3528688d956dca2cfbe1c95f46bd5a70998971 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 27 Feb 2023 11:28:12 -0700 Subject: [PATCH 19/44] Add scratch versions of entropy functions to base class --- singularity-eos/eos/eos_base.hpp | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index f91f62ba1b..1373f930a8 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -168,6 +168,13 @@ class EosBase { }); } 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, @@ -182,6 +189,13 @@ class EosBase { }); } 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, From b7ae939b13e6b2fc85c8bf332e7a709878cf44e9 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 27 Feb 2023 11:56:55 -0700 Subject: [PATCH 20/44] Ports of call update reverses condition to be more appropriate --- singularity-eos/eos/eos_ideal.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp index 7c0e6944cb..b7297bfd81 100644 --- a/singularity-eos/eos/eos_ideal.hpp +++ b/singularity-eos/eos/eos_ideal.hpp @@ -67,10 +67,10 @@ class IdealGas : public EosBase { 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_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 { From 8c65f1fc0af790af4c360709e9a5d62a20115fe7 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 27 Feb 2023 12:43:23 -0700 Subject: [PATCH 21/44] Fix a few typos when resolving merge conflicts --- singularity-eos/eos/eos_variant.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index 72ee4a88d7..d37bacc775 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -504,6 +504,8 @@ class Variant { [&rhos, &sies, &entropies, &num, &lambdas](const auto &eos) { return eos.EntropyFromDensityInternalEnergy(rhos, sies, entropies, num, lambdas); + }, + eos_); } template @@ -522,7 +524,7 @@ class Variant { RealIndexer &&entropies, Real *scratch, const int num, LambdaIndexer &&lambdas) const { return mpark::visit( - [&rhos, &sies, &pressures, &scratch, &num, &lambdas](const auto &eos) { + [&rhos, &sies, &entropies, &scratch, &num, &lambdas](const auto &eos) { return eos.EntropyFromDensityInternalEnergy(rhos, sies, entropies, scratch, num, lambdas); }, From d76d0f9fae5efd945b7ca0e664ea7104660c619f Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 27 Feb 2023 17:24:56 -0700 Subject: [PATCH 22/44] Fix extra template parameter --- singularity-eos/eos/eos_variant.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index d37bacc775..fa3ece4d96 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -420,7 +420,7 @@ class Variant { eos_); } - template + template inline void PressureFromDensityInternalEnergy(ConstRealIndexer &&rhos, ConstRealIndexer &&sies, RealIndexer &&pressures, Real *scratch, From 3d019ad293df023394801823776200b918db94dd Mon Sep 17 00:00:00 2001 From: Jeff Peterson Date: Mon, 27 Feb 2023 17:54:36 -0700 Subject: [PATCH 23/44] clang format --- singularity-eos/eos/eos_base.hpp | 19 ++++---- singularity-eos/eos/eos_ideal.hpp | 24 ++++------ singularity-eos/eos/eos_spiner.hpp | 16 +++---- singularity-eos/eos/eos_variant.hpp | 45 +++++++++---------- .../eos/modifiers/eos_unitsystem.hpp | 7 ++- singularity-eos/eos/modifiers/ramps_eos.hpp | 4 +- .../eos/modifiers/relativistic_eos.hpp | 4 +- singularity-eos/eos/modifiers/scaled_eos.hpp | 8 ++-- singularity-eos/eos/modifiers/shifted_eos.hpp | 4 +- 9 files changed, 61 insertions(+), 70 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 1373f930a8..f602839753 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -155,9 +155,9 @@ class EosBase { } template inline void EntropyFromDensityTemperature(ConstRealIndexer &&rhos, - ConstRealIndexer &&temperatures, - RealIndexer &&entropies, const int num, - LambdaIndexer &&lambdas) const { + 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)); @@ -170,8 +170,8 @@ class EosBase { template inline void EntropyFromDensityTemperature(ConstRealIndexer &&rhos, ConstRealIndexer &&temperatures, - RealIndexer &&entropies, Real * /*scratch*/, - const int num, LambdaIndexer &&lambdas) const { + RealIndexer &&entropies, Real * /*scratch*/, + const int num, LambdaIndexer &&lambdas) const { EntropyFromDensityTemperature(rhos, temperatures, entropies, num, lambdas); } template @@ -191,8 +191,8 @@ class EosBase { template inline void EntropyFromDensityInternalEnergy(ConstRealIndexer &&rhos, ConstRealIndexer &&sies, - RealIndexer &&entropies, Real * /*scratch*/, - const int num, LambdaIndexer &&lambdas) const { + RealIndexer &&entropies, Real * /*scratch*/, + const int num, LambdaIndexer &&lambdas) const { EntropyFromDensityInternalEnergy(rhos, sies, entropies, num, lambdas); } template @@ -386,9 +386,8 @@ class EosBase { Real RhoPmin(const Real temp) const { return 0.0; } // Default entropy behavior is to return an error - [[noreturn]] PORTABLE_FORCEINLINE_FUNCTION - void EntropyIsNotEnabled() const { - PORTABLE_ALWAYS_THROW_OR_ABORT(std::string("Entropy is not enabled for the '") + + [[noreturn]] PORTABLE_FORCEINLINE_FUNCTION void EntropyIsNotEnabled() const { + PORTABLE_ALWAYS_THROW_OR_ABORT(std::string("Entropy is not enabled for the '") + std::string(typeid(CRTP).name()) + std::string("' EOS")); } diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp index b7297bfd81..cc084c76c9 100644 --- a/singularity-eos/eos/eos_ideal.hpp +++ b/singularity-eos/eos/eos_ideal.hpp @@ -19,7 +19,6 @@ #include #include #include -#include // Ports-of-call #include @@ -42,24 +41,17 @@ 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)), _EntropyT0(_T0), - _EntropyRho0(_rho0) - { + _dvdt0(1. / (_rho0 * _T0)), _EntropyT0(_T0), _EntropyRho0(_rho0) { checkParams(); } - PORTABLE_INLINE_FUNCTION IdealGas(Real gm1, Real Cv, Real EntropyT0, - Real EntropyRho0) + 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) - { + _dvdt0(1. / (_rho0 * _T0)), _EntropyT0(EntropyT0), _EntropyRho0(EntropyRho0) { checkParams(); } - IdealGas GetOnDevice() { - return *this; - } + 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); @@ -69,8 +61,10 @@ class IdealGas : public EosBase { // 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_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 { @@ -87,7 +81,7 @@ class IdealGas : public EosBase { 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)); + _gm1 * _Cv * log(robust::ratio(_EntropyRho0, rho)); } PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda = nullptr) const { diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index c6678ccda8..0f02760d23 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -948,15 +948,15 @@ Real SpinerEOSDependsRhoT::PressureFromDensityInternalEnergy(const Real rho, PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::EntropyFromDensityTemperature(const Real rho, - const Real temperature, - Real *lambda) const { + const Real temperature, + Real *lambda) const { EntropyIsNotEnabled(); return 1.0; } PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::EntropyFromDensityInternalEnergy(const Real rho, - const Real sie, - Real *lambda) const { + const Real sie, + Real *lambda) const { EntropyIsNotEnabled(); return 1.0; } @@ -1679,15 +1679,15 @@ Real SpinerEOSDependsRhoSie::PressureFromDensityInternalEnergy(const Real rho, PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoSie::EntropyFromDensityTemperature(const Real rho, - const Real temperature, - Real *lambda) const { + const Real temperature, + Real *lambda) const { EntropyIsNotEnabled(); return 1.0; } PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoSie::EntropyFromDensityInternalEnergy(const Real rho, - const Real sie, - Real *lambda) const { + const Real sie, + Real *lambda) const { EntropyIsNotEnabled(); return 1.0; } diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index fa3ece4d96..02174cb5cb 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -446,29 +446,29 @@ class Variant { template inline void EntropyFromDensityTemperature(ConstRealIndexer &&rhos, ConstRealIndexer &&temperatures, - RealIndexer &&entropies, const int num) const { + 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 { + 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); + lambdas); }, eos_); } template - inline void - EntropyFromDensityTemperature(ConstRealIndexer &&rhos, ConstRealIndexer &&temperatures, - RealIndexer &&entropies, Real *scratch, - const int num) const { + 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); @@ -490,43 +490,42 @@ class Variant { template inline void EntropyFromDensityInternalEnergy(ConstRealIndexer &&rhos, ConstRealIndexer &&sies, - RealIndexer &&entropies, const int num) const { + 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 { + 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); + lambdas); }, eos_); } template inline void EntropyFromDensityInternalEnergy(ConstRealIndexer &&rhos, - ConstRealIndexer &&sies, - RealIndexer &&entropies, Real *scratch, - const int num) const { + 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); + 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 { + 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); + return eos.EntropyFromDensityInternalEnergy(rhos, sies, entropies, scratch, num, + lambdas); }, eos_); } diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index c1c2bad633..f6386f5dc0 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -85,8 +85,7 @@ class UnitSystem : public EosBase> { 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_entropy_unit_(1 / entropy_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? @@ -134,7 +133,7 @@ class UnitSystem : public EosBase> { } PORTABLE_FUNCTION Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, - Real *lambda = nullptr) const { + Real *lambda = nullptr) const { const Real S = t_.EntropyFromDensityInternalEnergy(rho * rho_unit_, sie * sie_unit_, lambda); return inv_entropy_unit_ * S; @@ -168,7 +167,7 @@ class UnitSystem : public EosBase> { return inv_press_unit_ * P; } Real EntropyFromDensityTemperature(const Real rho, const Real temp, - Real *lambda = nullptr) const { + Real *lambda = nullptr) const { const Real S = t_.EntropyFromDensityTemperature(rho * rho_unit_, temp * temp_unit_, lambda); return inv_entropy_unit_ * S; diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index ed668419a7..e3ab1e2663 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -144,7 +144,7 @@ class BilinearRampEOS : public EosBase> { } PORTABLE_FUNCTION Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, - Real *lambda = nullptr) const { + Real *lambda = nullptr) const { return t_.EntropyFromDensityInternalEnergy(rho, sie, lambda); } PORTABLE_FUNCTION @@ -175,7 +175,7 @@ class BilinearRampEOS : public EosBase> { } PORTABLE_FUNCTION Real EntropyFromDensityTemperature(const Real rho, const Real temperature, - Real *lambda = nullptr) const { + Real *lambda = nullptr) const { return t_.EntropyFromDensityTemperature(rho, temperature, lambda); } PORTABLE_FUNCTION diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp index a6d05f9242..80c7fbc347 100644 --- a/singularity-eos/eos/modifiers/relativistic_eos.hpp +++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp @@ -93,7 +93,7 @@ class RelativisticEOS : public EosBase> { } PORTABLE_FUNCTION Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, - Real *lambda = nullptr) const { + Real *lambda = nullptr) const { return t_.EntropyFromDensityInternalEnergy(rho, sie, lambda); } PORTABLE_FUNCTION @@ -121,7 +121,7 @@ class RelativisticEOS : public EosBase> { } PORTABLE_FUNCTION Real EntropyFromDensityTemperature(const Real rho, const Real temperature, - Real *lambda = nullptr) const { + Real *lambda = nullptr) const { return t_.EntropyFromDensityTemperature(rho, temperature, lambda); } PORTABLE_FUNCTION diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp index b0c854a15c..c71b600675 100644 --- a/singularity-eos/eos/modifiers/scaled_eos.hpp +++ b/singularity-eos/eos/modifiers/scaled_eos.hpp @@ -91,9 +91,9 @@ class ScaledEOS : public EosBase> { 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); + Real *lambda = nullptr) const { + return scale_ * + t_.EntropyFromDensityInternalEnergy(scale_ * rho, inv_scale_ * sie, lambda); } PORTABLE_FUNCTION Real SpecificHeatFromDensityInternalEnergy(const Real rho, const Real sie, @@ -120,7 +120,7 @@ class ScaledEOS : public EosBase> { } PORTABLE_FUNCTION Real EntropyFromDensityTemperature(const Real rho, const Real temperature, - Real *lambda = nullptr) const { + Real *lambda = nullptr) const { return scale_ * t_.EntropyFromDensityTemperature(scale_ * rho, temperature, lambda); } PORTABLE_FUNCTION diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index b165df1548..fde5182e0c 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -90,7 +90,7 @@ class ShiftedEOS : public EosBase> { } PORTABLE_FUNCTION Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, - Real *lambda = nullptr) const { + Real *lambda = nullptr) const { return t_.EntropyFromDensityInternalEnergy(rho, sie - shift_, lambda); } PORTABLE_FUNCTION @@ -115,7 +115,7 @@ class ShiftedEOS : public EosBase> { } PORTABLE_FUNCTION Real EntropyFromDensityTemperature(const Real rho, const Real temperature, - Real *lambda = nullptr) const { + Real *lambda = nullptr) const { return t_.EntropyFromDensityTemperature(rho, temperature, lambda); } PORTABLE_FUNCTION From bc96b66f9f33efcf7725a012c276b92febd21eb0 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 27 Feb 2023 18:55:50 -0700 Subject: [PATCH 24/44] Fix typo in ideal gas entropy lookup --- singularity-eos/eos/eos_ideal.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp index b7297bfd81..c8c0b2b6f7 100644 --- a/singularity-eos/eos/eos_ideal.hpp +++ b/singularity-eos/eos/eos_ideal.hpp @@ -19,7 +19,6 @@ #include #include #include -#include // Ports-of-call #include @@ -92,7 +91,7 @@ class IdealGas : public EosBase { 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, sie, lambda); + return EntropyFromDensityTemperature(rho, temp, lambda); } PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const { From 89f07d35ced1fb651e505c77cdbc3d400ec24b0d Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 27 Feb 2023 19:08:12 -0700 Subject: [PATCH 25/44] Add base class usings for entropy vector lookups --- singularity-eos/eos/eos_spiner.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index 0f02760d23..fa21f02c17 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; @@ -307,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; From 665174ed4c9489f02faa0cdf41795972a0cd0836 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 27 Feb 2023 19:11:46 -0700 Subject: [PATCH 26/44] Add base class usings to modifiers --- singularity-eos/eos/modifiers/eos_unitsystem.hpp | 2 ++ singularity-eos/eos/modifiers/relativistic_eos.hpp | 2 ++ singularity-eos/eos/modifiers/scaled_eos.hpp | 2 ++ singularity-eos/eos/modifiers/shifted_eos.hpp | 2 ++ 4 files changed, 8 insertions(+) diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index f6386f5dc0..3f59ceec34 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -57,6 +57,8 @@ class UnitSystem : public EosBase> { using EosBase>::PressureFromDensityInternalEnergy; using EosBase>::EntropyFromDensityTemperature; using EosBase>::EntropyFromDensityInternalEnergy; + using EosBase>::EntropyFromDensityTemperature; + using EosBase>::EntropyFromDensityInternalEnergy; using EosBase>::SpecificHeatFromDensityTemperature; using EosBase>::SpecificHeatFromDensityInternalEnergy; using EosBase>::BulkModulusFromDensityTemperature; diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp index 80c7fbc347..5ac26730a6 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; diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp index c71b600675..b2bbe78412 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; diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index fde5182e0c..fd303df543 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; From 408226ea56b75a351efe5c0ffcb1bcfac41c28ca Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 27 Feb 2023 19:12:37 -0700 Subject: [PATCH 27/44] Just kidding... this already had the usings --- singularity-eos/eos/modifiers/eos_unitsystem.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index 3f59ceec34..f6386f5dc0 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -57,8 +57,6 @@ class UnitSystem : public EosBase> { using EosBase>::PressureFromDensityInternalEnergy; using EosBase>::EntropyFromDensityTemperature; using EosBase>::EntropyFromDensityInternalEnergy; - using EosBase>::EntropyFromDensityTemperature; - using EosBase>::EntropyFromDensityInternalEnergy; using EosBase>::SpecificHeatFromDensityTemperature; using EosBase>::SpecificHeatFromDensityInternalEnergy; using EosBase>::BulkModulusFromDensityTemperature; From 818a495973a45dcbc5783e117a07d68c8d26c0df Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 27 Feb 2023 19:15:29 -0700 Subject: [PATCH 28/44] Add base class usings for entropy --- singularity-eos/eos/eos_stellar_collapse.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index da132df273..697aabd9cc 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; From f5cafe9993ca4d536a05c63dee3901d480d44575 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 27 Feb 2023 19:25:56 -0700 Subject: [PATCH 29/44] Needed to add a S(rho,e) method for stellar collapse because there was just S(rho,T) --- singularity-eos/eos/eos_stellar_collapse.hpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 697aabd9cc..0ad67f24f9 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -109,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 @@ -505,6 +508,16 @@ 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, From 240198c2a39d075f4c77b4b9b41f69bd56fd731e Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 27 Feb 2023 19:31:39 -0700 Subject: [PATCH 30/44] Add default error messages for JWL and Davis when accessing entropy (for now at least) --- singularity-eos/eos/eos_davis.hpp | 10 ++++++++++ singularity-eos/eos/eos_jwl.hpp | 14 ++++++++++++++ 2 files changed, 24 insertions(+) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 1c151ade58..58aff6520d 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -162,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(); + return 1.0; + } + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + EntropyIsNotEnabled(); + 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_jwl.hpp b/singularity-eos/eos/eos_jwl.hpp index 977608b145..1bff918f63 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,11 @@ 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(); + 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 +155,11 @@ 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(); + return 1.0; +} PORTABLE_INLINE_FUNCTION Real JWL::SpecificHeatFromDensityTemperature( const Real rho, const Real temp, Real *lambda) const { return SpecificHeatFromDensityInternalEnergy( From 1a1df04d13749df730cb953e985780ffe70304b3 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 27 Feb 2023 19:45:33 -0700 Subject: [PATCH 31/44] Add basic test for ideal gas entropy --- test/test_eos_unit.cpp | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index f916140292..a5ed5c3e24 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -517,15 +517,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 +557,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 +595,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); From a5ea8ec7a8aed9433027747da6d4d358755522bb Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 28 Feb 2023 09:52:36 -0700 Subject: [PATCH 32/44] Add scalar test for ideal gas entropy at key values --- test/test_eos_unit.cpp | 75 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index a5ed5c3e24..02f4f28b38 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -467,6 +467,44 @@ SCENARIO("EOS Unit System", "[EOSBuilder][UnitSystem][IdealGas]") { } } +SCENARIO("Ideal 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)") { + constexpr 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)") { + constexpr 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") { @@ -687,15 +725,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); @@ -704,12 +765,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(); @@ -739,6 +802,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); From 5e5ca0b4d76c04c6ea26ad7e8e957f320c0b80b6 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 28 Feb 2023 10:30:11 -0700 Subject: [PATCH 33/44] Add test of the entropy error in the Gruneisen EOS --- test/test_eos_gruneisen.cpp | 28 ++++++++++++++++++++++++++++ test/test_eos_unit.cpp | 3 +-- 2 files changed, 29 insertions(+), 2 deletions(-) diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index caaf6e1947..f07ba9e63d 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -32,6 +32,34 @@ PORTABLE_INLINE_FUNCTION Real QuadFormulaMinus(Real a, Real b, Real c) { return (-b - std::sqrt(b * b - 4 * a * c)) / (2 * a); } +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); + } + } +} + 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 02f4f28b38..5b95b0beb6 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -467,8 +467,7 @@ SCENARIO("EOS Unit System", "[EOSBuilder][UnitSystem][IdealGas]") { } } -SCENARIO("Ideal entropy", "[IdealGas][Entropy]") { - +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; From 743d1dc6849c719a847909b1626766873ec8a6ee Mon Sep 17 00:00:00 2001 From: Jeff Peterson Date: Tue, 28 Feb 2023 10:33:17 -0700 Subject: [PATCH 34/44] clang format --- singularity-eos/eos/eos_jwl.hpp | 10 ++++--- singularity-eos/eos/eos_stellar_collapse.hpp | 5 ++-- test/test_eos_unit.cpp | 30 +++++++++++--------- 3 files changed, 24 insertions(+), 21 deletions(-) diff --git a/singularity-eos/eos/eos_jwl.hpp b/singularity-eos/eos/eos_jwl.hpp index 1bff918f63..566d242b7c 100644 --- a/singularity-eos/eos/eos_jwl.hpp +++ b/singularity-eos/eos/eos_jwl.hpp @@ -122,8 +122,9 @@ 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 { +PORTABLE_INLINE_FUNCTION Real JWL::EntropyFromDensityInternalEnergy(const Real rho, + const Real sie, + Real *lambda) const { EntropyIsNotEnabled(); return 1.0; } @@ -155,8 +156,9 @@ 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 { +PORTABLE_INLINE_FUNCTION Real JWL::EntropyFromDensityTemperature(const Real rho, + const Real temp, + Real *lambda) const { EntropyIsNotEnabled(); return 1.0; } diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 0ad67f24f9..eeb2d3bbf7 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -110,7 +110,7 @@ class StellarCollapse : public EosBase { Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie, - Real *lambda = nullptr) const; + Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(const Real rho, const Real temperature, Real *lambda = nullptr) const; @@ -509,8 +509,7 @@ Real StellarCollapse::EntropyFromDensityTemperature(const Real rho, } PORTABLE_INLINE_FUNCTION -Real StellarCollapse::EntropyFromDensityInternalEnergy(const Real rho, - const Real sie, +Real StellarCollapse::EntropyFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda) const { Real lRho, lT, Ye; getLogsFromRhoSie_(rho, sie, lambda, lRho, lT, Ye); diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index 5b95b0beb6..30010ec83f 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -481,9 +481,10 @@ SCENARIO("Ideal gas entropy", "[IdealGas][Entropy]") { 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") { + 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 + constexpr Real rho = 0.1; // rho**3 = EntropyRho0 THEN("The entropy should be 2. / 3. * gm1 * Cv * log(EntropyRho0)") { constexpr Real entropy_true = 2. / 3. * gm1 * Cv * log(EntropyRho0); auto entropy = eos.EntropyFromDensityTemperature(rho, T); @@ -491,8 +492,9 @@ SCENARIO("Ideal gas entropy", "[IdealGas][Entropy]") { 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 + 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)") { constexpr Real entropy_true = -1. / 2. * Cv * log(EntropyT0); @@ -562,14 +564,14 @@ SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") { std::array entropy_true; Real *v_entropy_true = entropy_true.data(); #endif - constexpr Real P0 = 1e6; // microbar - constexpr Real T0 = 293; // K + 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]); - }); + 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); @@ -732,14 +734,14 @@ SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") { std::array entropy_true; Real *v_entropy_true = entropy_true.data(); #endif - constexpr Real P0 = 1e6; // microbar - constexpr Real T0 = 293; // K + 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]); - }); + 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); From 4c3d2aee3e11761e11bad483cb6f4028ddc4635d Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 28 Feb 2023 12:42:48 -0700 Subject: [PATCH 35/44] Update ports-of-call in the package.py to update the version in the gitlab CI so that it passes --- spack-repo/packages/singularity-eos/package.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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") From 8382c4b735ceed3eae2ef8d0ba3255cd2fdbbd4c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 28 Feb 2023 19:15:53 -0700 Subject: [PATCH 36/44] Change the EntropyIsNotEnabled() function to use char and remove typeid introspection to allow to work on device --- singularity-eos/eos/eos_base.hpp | 12 ++++++++---- singularity-eos/eos/eos_davis.hpp | 8 ++++---- singularity-eos/eos/eos_eospac.hpp | 2 +- singularity-eos/eos/eos_gruneisen.hpp | 2 +- singularity-eos/eos/eos_jwl.hpp | 4 ++-- singularity-eos/eos/eos_spiner.hpp | 8 ++++---- 6 files changed, 20 insertions(+), 16 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index f602839753..09e202045a 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -16,6 +16,7 @@ #define _SINGULARITY_EOS_EOS_EOS_BASE_ #include +#include #include #include @@ -386,10 +387,13 @@ class EosBase { Real RhoPmin(const Real temp) const { return 0.0; } // Default entropy behavior is to return an error - [[noreturn]] PORTABLE_FORCEINLINE_FUNCTION void EntropyIsNotEnabled() const { - PORTABLE_ALWAYS_THROW_OR_ABORT(std::string("Entropy is not enabled for the '") + - std::string(typeid(CRTP).name()) + - std::string("' EOS")); + 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 '"; + strcat(msg, eosname); + strcat(msg, "' EOS"); + PORTABLE_ALWAYS_THROW_OR_ABORT(msg); } }; } // namespace eos_base diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 58aff6520d..82e54e9ff5 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -62,12 +62,12 @@ class DavisReactants : public EosBase { } PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const { - EntropyIsNotEnabled(); + EntropyIsNotEnabled("DavisReactants"); return 1.0; } PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda = nullptr) const { - EntropyIsNotEnabled(); + EntropyIsNotEnabled("DavisReactants"); return 1.0; } PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( @@ -164,12 +164,12 @@ class DavisProducts : public EosBase { } PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const { - EntropyIsNotEnabled(); + EntropyIsNotEnabled("DavisProducts"); return 1.0; } PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda = nullptr) const { - EntropyIsNotEnabled(); + EntropyIsNotEnabled("DavisProducts"); return 1.0; } PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp index c828165ed7..67bdda5a89 100644 --- a/singularity-eos/eos/eos_eospac.hpp +++ b/singularity-eos/eos/eos_eospac.hpp @@ -819,7 +819,7 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::PressureFromDensityTemperature(const Real PORTABLE_INLINE_FUNCTION Real EOSPAC::EntropyFromDensityTemperature( const Real rho, const Real temperature, Real *lambda) const { - EntropyIsNotEnabled(); + EntropyIsNotEnabled("EOSPAC"); return 1.0; } diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index 385cce9051..d1fcc3622a 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -317,7 +317,7 @@ PORTABLE_INLINE_FUNCTION Real Gruneisen::PressureFromDensityInternalEnergy( 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(); + EntropyIsNotEnabled("Gruneisen"); return 1.0; } PORTABLE_INLINE_FUNCTION Real Gruneisen::BulkModulusFromDensityInternalEnergy( diff --git a/singularity-eos/eos/eos_jwl.hpp b/singularity-eos/eos/eos_jwl.hpp index 566d242b7c..201d33b00a 100644 --- a/singularity-eos/eos/eos_jwl.hpp +++ b/singularity-eos/eos/eos_jwl.hpp @@ -125,7 +125,7 @@ PORTABLE_INLINE_FUNCTION Real JWL::PressureFromDensityInternalEnergy(const Real PORTABLE_INLINE_FUNCTION Real JWL::EntropyFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda) const { - EntropyIsNotEnabled(); + EntropyIsNotEnabled("JWL"); return 1.0; } PORTABLE_INLINE_FUNCTION Real JWL::TemperatureFromDensityInternalEnergy( @@ -159,7 +159,7 @@ PORTABLE_INLINE_FUNCTION Real JWL::PressureFromDensityTemperature(const Real rho PORTABLE_INLINE_FUNCTION Real JWL::EntropyFromDensityTemperature(const Real rho, const Real temp, Real *lambda) const { - EntropyIsNotEnabled(); + EntropyIsNotEnabled("JWL"); return 1.0; } PORTABLE_INLINE_FUNCTION Real JWL::SpecificHeatFromDensityTemperature( diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index fa21f02c17..fdb9c672a8 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -954,14 +954,14 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::EntropyFromDensityTemperature(const Real rho, const Real temperature, Real *lambda) const { - EntropyIsNotEnabled(); + EntropyIsNotEnabled("SpinerEOSDependsRhoT"); return 1.0; } PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::EntropyFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda) const { - EntropyIsNotEnabled(); + EntropyIsNotEnabled("SpinerEOSDependsRhoT"); return 1.0; } @@ -1685,14 +1685,14 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoSie::EntropyFromDensityTemperature(const Real rho, const Real temperature, Real *lambda) const { - EntropyIsNotEnabled(); + EntropyIsNotEnabled("SpinerEOSDependsRhoSie"); return 1.0; } PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoSie::EntropyFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda) const { - EntropyIsNotEnabled(); + EntropyIsNotEnabled("SpinerEOSDependsRhoSie"); return 1.0; } From 7c2452b408513f4e6215c63c793a860b5c4842d1 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 28 Feb 2023 20:06:49 -0700 Subject: [PATCH 37/44] Fix some warnings --- doc/sphinx/src/sphinx-doc.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) 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 From 2f985bb3445010b795ad960563580faaa8d055bc Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 28 Feb 2023 20:07:14 -0700 Subject: [PATCH 38/44] Add some footnotes and update the Davis EOS --- doc/sphinx/src/models.rst | 178 +++++++++++++++++++++++--------------- 1 file changed, 110 insertions(+), 68 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index fd7bd206e5..6d1b5b7da8 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -73,13 +73,13 @@ 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. -Completing an Incomplete Equation of State and Consequences for Thermodynamic Consistency -````````````````````````````````````````````````````````````````````````````````````````` +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 +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. @@ -95,13 +95,13 @@ variables*, i.e. F =& F(\rho, T) \\ G =& G(P, T) \\ -where all potentials are specific, :math:`e` is again the internal energy, +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, the Mie Gruneisen -form tends to dominate most analytic EOS and it becomes imperitive to extend the -incomplete form to a complete equation of state. +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 @@ -165,42 +165,52 @@ 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, 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)` | Density | Pressure, Temperature | -+--------------------------+----------------------+ | -| :math:`e(P, T)` | 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 | -+--------------------------+----------------------+--------------------------+ ++---------------------------+----------------------+--------------------------+ +| 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 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. Thus + for a multivalued pressure, the returned density will most often be that + closest to the STP density. 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 and entopry @@ -673,7 +683,7 @@ 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 @@ -702,7 +712,7 @@ Davis Reactants EOS .. warning:: Entropy is not yet available for this EOS -The `Davis reactants EOS ` uses an isentrope passing through a +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 @@ -713,37 +723,77 @@ capacity varies linearly with entropy such that where subscript :math:`0` refers to the reference state and :math:`\alpha` is a dimensionless constant specified by the user. -Using the fact that the heat capacity is defined by - -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 -The Gruneisen parameter takes on a linear form such that +.. math:: + + \Gamma_S(\rho) = \frac{\rho}{T} + \left(\frac{\partial T}{\partial \rho}\right)_S, + +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` is a dimensionless parameter and :math:`y = 1 - \rho0/\rho`. +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} @@ -751,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` From 411cec64c2b4537203248f629ab7472976b95d8b Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 28 Feb 2023 20:09:32 -0700 Subject: [PATCH 39/44] Remove constexpr to satisfy power9 compiler --- test/test_eos_unit.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index 30010ec83f..ba46084755 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -486,7 +486,7 @@ SCENARIO("Ideal gas entropy", "[IdealGas][Entropy]") { constexpr Real T = EntropyT0; constexpr Real rho = 0.1; // rho**3 = EntropyRho0 THEN("The entropy should be 2. / 3. * gm1 * Cv * log(EntropyRho0)") { - constexpr Real entropy_true = 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)); @@ -497,7 +497,7 @@ SCENARIO("Ideal gas entropy", "[IdealGas][Entropy]") { constexpr Real T = 10; // T**2 = EntropyT0 constexpr Real rho = EntropyRho0; THEN("The entropy should be -1. / 2. * Cv * log(EntropyT0)") { - constexpr Real entropy_true = -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)); From 3fc6125e71d620c8f8183de0702563506bea4d2d Mon Sep 17 00:00:00 2001 From: Jeff Peterson Date: Tue, 28 Feb 2023 20:14:22 -0700 Subject: [PATCH 40/44] clang format --- singularity-eos/eos/eos_base.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 09e202045a..fded3e9c30 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -15,8 +15,8 @@ #ifndef _SINGULARITY_EOS_EOS_EOS_BASE_ #define _SINGULARITY_EOS_EOS_EOS_BASE_ -#include #include +#include #include #include @@ -388,7 +388,7 @@ class EosBase { // Default entropy behavior is to return an error PORTABLE_FORCEINLINE_FUNCTION - void EntropyIsNotEnabled(const char* eosname) const { + 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 '"; strcat(msg, eosname); From e220b751dcf59dc69910213663a0ca3d644f1b8e Mon Sep 17 00:00:00 2001 From: Jeff Peterson <83598606+jhp-lanl@users.noreply.github.com> Date: Wed, 1 Mar 2023 12:56:07 -0700 Subject: [PATCH 41/44] Add `std` namespace to srcat calls Co-authored-by: Jonah Miller --- singularity-eos/eos/eos_base.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index fded3e9c30..b413245212 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -391,8 +391,8 @@ class EosBase { 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 '"; - strcat(msg, eosname); - strcat(msg, "' EOS"); + std::strcat(msg, eosname); + std::strcat(msg, "' EOS"); PORTABLE_ALWAYS_THROW_OR_ABORT(msg); } }; From 5b894e05c8fa834f55ea6e0cf9c4d7a663b38d89 Mon Sep 17 00:00:00 2001 From: Jeff Peterson <83598606+jhp-lanl@users.noreply.github.com> Date: Wed, 1 Mar 2023 16:37:41 -0700 Subject: [PATCH 42/44] Update language on P-T inversion to note spiner exception --- doc/sphinx/src/models.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 6d1b5b7da8..a82654e17b 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -205,12 +205,12 @@ in the :doc:`EOS API `. .. [#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 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. Thus - for a multivalued pressure, the returned density will most often be that - closest to the STP density. + 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 and entopry From 8c870cebe43cc362c542932d5be48f55ba59625e Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 1 Mar 2023 16:56:50 -0700 Subject: [PATCH 43/44] Only run exception test when portability strategy is none --- test/test_eos_gruneisen.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index f07ba9e63d..a62dc6f194 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -32,6 +32,8 @@ 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 @@ -59,6 +61,7 @@ SCENARIO("Gruneisen EOS entropy is disabled", "[GruneisenEOS][Entropy]") { } } } +#endif SCENARIO("Gruneisen EOS", "[VectorEOS][GruneisenEOS]") { GIVEN("Parameters for a Gruneisen EOS") { From 4655c062baaa8c1ed71dc8b711ed8bf2ed064787 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 16 Mar 2023 11:24:38 -0600 Subject: [PATCH 44/44] Cite Ann's reference --- doc/sphinx/src/models.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index a82654e17b..f9b15d3c0f 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -84,9 +84,9 @@ 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. -In general, an equation of state can be considered "complete" if it is -constructed from one of the thermodynamic potentials using their *natural -variables*, i.e. +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::