Added documentation for the simple MACAW some more comments.

Author: bgclayton-lanl

doc/sphinx/src/models.rst

@@ -21,6 +21,8 @@
 
 .. _PowerMG: https://www.osti.gov/biblio/1762624
 
+.. _SimpleMACAW: https://www.osti.gov/biblio/2479474
+
 
 EOS Models
 ===========
@@ -1320,6 +1322,58 @@ This constructor also optionally accepts `MeanAtomicProperties` for
 the atomic mass and number as a final optional parameter.
 
 
+.. _MACAW EOS: https://pubs.aip.org/aip/jap/article/134/12/125102/2912256
+
+Simple MACAW
+````````````
+The `Simple MACAW EOS <_SimpleMACAW>`_ is a simplified version of the `MACAW EOS <MACAW EOS>`_
+and is thermodynamically consistent.  It is constructed from a the Helmholtz
+free energy using a Murnaghan cold curve and a simplified thermal component
+from the MACAW EOS.
+
+Fundamentally, the equation of state can be written in Mie-Gruneisen form (with constant Gruneisen parameter) as:
+
+.. math::
+
+    P(v, e) = P_{\text{cold}}(v) + \Gamma_c \rho (e - e_{\text{cold}}(v))
+
+where the cold curves are given by:
+
+.. math::
+
+    e_{\text{cold}}(v) = A v_0 \Big[ \Big( \frac{v}{v_0} \Big)^{-B} + \Big( \frac{v}{v_0} \Big) B - (B+1) \Big]
+
+and 
+
+.. math::
+
+    p_{\text{cold}}(v) := -e'_{\text{cold}}(v) = AB \Big( \Big( \frac{v}{v_0} \Big)^{-(B+1)} - 1 \Big)
+
+The specific heat capacity at constant volume for this model is given by,
+
+.. math::
+
+    C_v(v, \tau) = C^\infty_v \frac{\tau^2 + 2\tau}{(\tau + 1)^2} \quad \text{where } \tau = \frac{T}{\theta(v)} \quad \text{ and } \quad \theta(v) := T_0 \Big( \frac{v}{v_c} \Big)^{-\Gamma_c}
+
+Note it obeys the expected physical behavior; that, :math:`\lim_{T\to 0^+} C_v(v,\tau(v,T)) = 0` and
+:math:`\lim_{T\to\infty} C_v(v,\tau(v,T) = C^\infty_v < \infty` (Dulong-Petit law).
+
+The constructor for the Simple MACAW EOS is
+
+.. code-block:: cpp
+
+  SimpleMACAW(const Real A, const Real B, const Real Cvinf, const Real v0,
+              const Real T0, const Real Gc)
+
+where ``A`` is :math:`A`, ``B`` is :math:`B`, ``Cvinf`` is :math:`C^\infty_v`,
+``v0`` is :math:`v_0`, ``T0`` is :math:`T_0`, ``Gc`` is :math:`\Gamma_c`.
+
+In order to maintain thermodynamic consistency, a sufficient set of constraints
+is given by :math:`A > 0`, :math:`B > 0`, :math:`C^\infty_v > 0`, :math:`v_0 >
+0`, :math:`T_0 > 0`, and :math:`\Gamma_c > 0`. If one wants to maintain
+thermodynamic consistency, we also require :math:`\Gamma_c \in (0,1)`. 
+
+
 JWL EOS
 ``````````
 

singularity-eos/eos/eos_simple_macaw.hpp

@@ -33,8 +33,12 @@ namespace singularity {
 
 using namespace eos_base;
 
-/* The details of this equation of state can be found here:
- * https://www.osti.gov/biblio/2479474 */
+/* Most information regarding this equation of state can be found here:
+ * https://www.osti.gov/biblio/2479474 
+ * 
+ * Note that all equations are given as functions of the density ($\rho$) rather than functions
+ * of the specific volume ($v$). The reason for this is it simplifies the computational complexity.
+ */
 class SimpleMACAW : public EosBase<SimpleMACAW> {
  public:
   SimpleMACAW() = default;
@@ -55,10 +59,10 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
     return robust::ratio((Delta_e + std::sqrt(discriminant)), 2.0 * _Cvinf);
   }
   PORTABLE_INLINE_FUNCTION void CheckParams() const {
-    PORTABLE_ALWAYS_REQUIRE(_A > 0, "Parameter, 'A', must be positive");
-    PORTABLE_ALWAYS_REQUIRE(_B > 0, "Parameter, 'B', must be positive");
+    PORTABLE_ALWAYS_REQUIRE(_A >= 0, "Parameter, 'A', must be non-negative");
+    PORTABLE_ALWAYS_REQUIRE(_B >= 0, "Parameter, 'B', must be non-negative");
     PORTABLE_ALWAYS_REQUIRE(_v0 > 0, "Reference specific volume, 'v0', must be positive");
-    PORTABLE_ALWAYS_REQUIRE(_T0 > 0, "Reference temperature, 'T0', must be positive");
+    PORTABLE_ALWAYS_REQUIRE(_T0 >= 0, "Reference temperature, 'T0', must be non-negative");
     PORTABLE_ALWAYS_REQUIRE(_Cvinf > 0, "Specific heat capacity (at constant volume), `Cvinf`, must be positive");
     PORTABLE_ALWAYS_REQUIRE(_Gc > 0 && _Gc < 1, "Gruneisen parameter, 'Gc', must be in the interval (0,1)");
     _AZbar.CheckParams();
@@ -67,18 +71,21 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
   PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature(
       const Real rho, const Real temperature,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+      /* Equation (16) */
     return rho * (SieColdCurve(rho) + SieThermalPortion(rho, temperature));
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature(
       const Real rho, const Real temperature,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+      /* Equation (15) */
     return PressureColdCurve(rho) + PressureThermalPortion(rho, temperature);
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy(
       const Real rho, const Real sie,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+      /* Equation (19) */
     return PressureColdCurve(rho) + _Gc * rho * (sie - SieColdCurve(rho));
   }
   template <typename Indexer_t = Real *>
@@ -93,6 +100,7 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
   PORTABLE_INLINE_FUNCTION Real
   EntropyFromDensityTemperature(const Real rho, const Real temperature,
                                 Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    /* Equation (8) */
     const Real tau = robust::ratio(temperature, TemperatureScale(rho));
     return _Cvinf * (robust::ratio(tau, 1.0 + tau) + std::log(1.0 + tau));
   }
@@ -138,6 +146,7 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real TemperatureScale(
       const Real rho, Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    /* Equation (13) */
     return _T0 * std::pow(rho * _v0, _Gc);
   }
 
@@ -146,6 +155,7 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
   PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(
       const Real rho, const Real temperature,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    /* Equation (6) */
     const Real tau = robust::ratio(temperature, TemperatureScale(rho));
     const Real numerator = math_utils::pow<2>(tau) + 2.0 * tau;
     const Real denominator = math_utils::pow<2>(tau + 1.0);
@@ -170,6 +180,7 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
   PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy(
       const Real rho, const Real sie,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    /* Equation (21) (multiplied by $\rho$ to get bulk mod) */
     const Real term1 = _A * _B * (_B + 1.0) * std::pow(rho * _v0, _B + 1.0);
     const Real term2 = _Gc * (_Gc + 1.0) * rho * (sie - SieColdCurve(rho));
     return term1 + term2;
@@ -179,7 +190,8 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
   PORTABLE_INLINE_FUNCTION Real IsothermalBulkModulusFromDensityTemperature(
       const Real rho, const Real temperature,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
-    /* As mentioned in the paper, from the constraints on the parameters,
+    /* Equation (22) (multiplied by $\rho$ to get bulk mod)
+     * As mentioned in the paper, from the constraints on the parameters,
      * the isothermal bulk modulus is guaranteed to be positive. */
     const Real term1 = _A * _B * (_B + 1.0) * std::pow(rho * _v0, _B + 1.0);
     const Real numerator = rho * _T0 * (1.0 - _Gc) * std::pow(rho * _v0, 2.0 * _Gc) + temperature;