I think I got all the corrections done.

Author: bgclayton-lanl

doc/sphinx/src/models.rst

@@ -1367,11 +1367,11 @@ The constructor for the Simple MACAW EOS is
 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
+In order to maintain thermodynamic stability, 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)`. 
-
+0`, :math:`T_0 > 0`, and :math:`\Gamma_c \in (0,1]`. One can still select
+:math:`\Gamma_c > 1`, just note that the isothermal bulk modulus can be
+negative (the isentropic bulk modulus will still be positive though).
 
 JWL EOS
 ``````````

singularity-eos/eos/eos_simple_macaw.hpp

@@ -45,7 +45,7 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
   PORTABLE_INLINE_FUNCTION
   SimpleMACAW(Real A, Real B, Real Cvinf, Real v0, Real T0, Real Gc,
            const MeanAtomicProperties &AZbar = MeanAtomicProperties())
-      : _A(A), _B(B), _Cvinf(Cvinf), _v0(v0), _T0(T0), _Gc(Gc), _AZbar(AZbar) {
+      : A_(A), B_(B), Cvinf_(Cvinf), v0_(v0), T0_(T0), Gc_(Gc), _AZbar(AZbar) {
     CheckParams();
   }
   SimpleMACAW GetOnDevice() { return *this; }
@@ -55,45 +55,45 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
     /* Equation (18) */
     const Real Delta_e = sie - SieColdCurve(rho);
-    const Real discriminant = Delta_e * (Delta_e + 4.0 * _Cvinf * _T0 * std::pow(rho * _v0, _Gc));
-    return robust::ratio((Delta_e + std::sqrt(discriminant)), 2.0 * _Cvinf);
+    const Real discriminant = Delta_e * (Delta_e + 4.0 * Cvinf_ * T0_ * std::pow(rho * v0_, Gc_));
+    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 non-negative");
-    PORTABLE_ALWAYS_REQUIRE(_B >= 0, "Parameter, 'B', must be non-negative");
-    PORTABLE_ALWAYS_REQUIRE(_Cvinf > 0, "Specific heat capacity (at constant volume), `Cvinf`, must be positive");
-    PORTABLE_ALWAYS_REQUIRE(_v0 > 0, "Reference specific volume, 'v0', must be positive");
-    PORTABLE_ALWAYS_REQUIRE(_T0 >= 0, "Reference temperature, 'T0', must be non-negative");
-    PORTABLE_ALWAYS_REQUIRE(_Gc > 0, "Gruneisen parameter, 'Gc', must be positive");
-    if (_Gc > 1) {PORTABLE_WARN("Warning: Gruneisen coefficient greater than 1. Thermodynamic stability may be violated.");}
+    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(Cvinf_ > 0, "Specific heat capacity (at constant volume), `Cvinf`, must be positive");
+    PORTABLE_ALWAYS_REQUIRE(v0_ > 0, "Reference specific volume, 'v0', must be positive");
+    PORTABLE_ALWAYS_REQUIRE(T0_ >= 0, "Reference temperature, 'T0', must be non-negative");
+    PORTABLE_ALWAYS_REQUIRE(Gc_ > 0, "Gruneisen parameter, 'Gc', must be positive");
+    if (Gc_ > 1) {PORTABLE_WARN("Warning: Gruneisen coefficient greater than 1. Thermodynamic stability may be violated.");}
     _AZbar.CheckParams();
   }
   template <typename Indexer_t = Real *>
   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));
+    return 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);
+    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));
+    return PressureColdCurve(rho) + Gc_ * rho * (sie - SieColdCurve(rho));
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityPressure(
       const Real rho, const Real P,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
-    return SieColdCurve(rho) + robust::ratio(P - PressureColdCurve(rho), _Gc * rho);
+    return SieColdCurve(rho) + robust::ratio(P - PressureColdCurve(rho), Gc_ * rho);
   }
 
   // Entropy member functions
@@ -102,8 +102,8 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
   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));
+    const Real tau = robust::ratio(temperature, _TemperatureScale(rho));
+    return Cvinf_ * (robust::ratio(tau, 1.0 + tau) + std::log(1.0 + tau));
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy(
@@ -120,35 +120,14 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real SieColdCurve(
       const Real rho, Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
-    const Real ratio = rho * _v0;
-    return _A * _v0 * (std::pow(ratio, _B) + robust::ratio(_B, ratio) - (_B + 1.));
-  }
-  template <typename Indexer_t = Real *>
-  PORTABLE_INLINE_FUNCTION Real SieThermalPortion(
-      const Real rho, const Real T,
-      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
-    const Real ratio = rho * _v0;
-    return _Cvinf * T * (1.0 - robust::ratio(_T0, _T0 + T * std::pow(ratio, -_Gc)));
+    const Real ratio = rho * v0_;
+    return A_ * v0_ * (std::pow(ratio, B_) + robust::ratio(B_, ratio) - (B_ + 1.));
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real PressureColdCurve(
       const Real rho, Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
-    const Real ratio = rho * _v0;
-    return _A * _B * (std::pow(ratio, _B + 1.0) - 1.0);
-  }
-  template <typename Indexer_t = Real *>
-  PORTABLE_INLINE_FUNCTION Real PressureThermalPortion(
-      const Real rho, const Real T,
-      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
-      const Real numerator = rho * _Cvinf * _Gc * math_utils::pow<2>(T);
-      const Real denominator = _T0 * std::pow(rho * _v0, _Gc) + T;
-    return robust::ratio(numerator, denominator);
-  }
-  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);
+    const Real ratio = rho * v0_;
+    return A_ * B_ * (std::pow(ratio, B_ + 1.0) - 1.0);
   }
 
   // Specific heat capacity at constant volume
@@ -157,10 +136,10 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
       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 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);
-    return _Cvinf * robust::ratio(numerator, denominator);
+    return Cvinf_ * robust::ratio(numerator, denominator);
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy(
@@ -182,8 +161,8 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
       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));
+    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;
   }
   // Isothermal bulk modulus
@@ -194,10 +173,10 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
     /* 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;
-    const Real denominator = _T0 * std::pow(rho * _v0, _Gc) + temperature;
-    return term1 + _Cvinf * _Gc * temperature * temperature * robust::ratio(numerator, denominator * denominator);
+    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;
+    const Real denominator = T0_ * std::pow(rho * v0_, Gc_) + temperature;
+    return term1 + Cvinf_ * Gc_ * temperature * temperature * robust::ratio(numerator, denominator * denominator);
   }
   // Specific heat capacity at constant pressure
    template <typename Indexer_t = Real *>
@@ -215,28 +194,26 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
   PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature(
       const Real rho, const Real temperature,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
-    return _Gc;
+    return Gc_;
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy(
       const Real rho, const Real sie,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
-    return _Gc;
+    return Gc_;
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod,
                   const unsigned long output, Indexer_t &&lambda) const {
   if (output & thermalqs::density && output & thermalqs::specific_internal_energy) {
-    if (output & thermalqs::pressure || output & thermalqs::temperature) {
-      UNDEFINED_ERROR;
-    }
+    PORTABLE_ALWAYS_REQUIRE(!(output & thermalqs::pressure || output & thermalqs::temperature), 
+                            "Cannot request more than two thermodynamic variables (rho, T, e, P)");
     DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie);
   }
   if (output & thermalqs::pressure && output & thermalqs::specific_internal_energy) {
-    if (output & thermalqs::density || output & thermalqs::temperature) {
-      UNDEFINED_ERROR;
-    }
+    PORTABLE_ALWAYS_REQUIRE(!(output & thermalqs::density || output & thermalqs::temperature), 
+                            "Cannot request more than two thermodynamic variables (rho, T, e, P)");
     sie = InternalEnergyFromDensityTemperature(rho, temp, lambda);
   }
   if (output & thermalqs::temperature && output & thermalqs::specific_internal_energy) {
@@ -256,8 +233,8 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
                          Real &bmod, Real &dpde, Real &dvdt,
                          Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
     // v_0 and T_0 are the only user selected reference state parameters
-    rho = 1.0 / _v0;
-    temp = _T0;
+    rho = 1.0 / v0_;
+    temp = T0_;
     press = PressureFromDensityTemperature(rho, temp);
     sie = InternalEnergyFromDensityTemperature(rho, temp);
 
@@ -274,7 +251,7 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
   PORTABLE_INLINE_FUNCTION void PrintParams() const {
     printf("Simple MACAW Parameters:\nA = %g\nB    = %g\nCvinf  = %g\nv0  = %g\nT0  = %g\nGc     = "
            "%g\n",
-           _A, _B, _Cvinf, _v0, _T0, _Gc);
+           A_, B_, Cvinf_, v0_, T0_, Gc_);
     _AZbar.PrintParams();
   }
   template <typename Indexer_t>
@@ -286,19 +263,19 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
 
     // Setup lambda function for rho
     auto f = [=](const Real x /* density */){
-      const Real term1 = _A * _B * (std::pow(_v0 * x, _B + 1.0) - 1.0) - press;
-      const Real term2 = _T0 * std::pow(_v0 * x, _Gc) + temp;
-      const Real term3 = _Cvinf * _Gc * temp * temp * x;
+      const Real term1 = A_ * B_ * (std::pow(v0_ * x, B_ + 1.0) - 1.0) - press;
+      const Real term2 = T0_ * std::pow(v0_ * x, Gc_) + temp;
+      const Real term3 = Cvinf_ * Gc_ * temp * temp * x;
       return term1 * term2 + term3;
     };
 
     const RootFinding1D::RootCounts root_info;
 
     // Finding the density on pressure cold curve is a guaranteed upper bound
-    const Real rho_high = std::pow(press / (_A * _B) + 1.0, 1.0 / (_B + 1.0)) / _v0;
+    const Real rho_high = std::pow(press / (A_ * B_) + 1.0, 1.0 / (B_ + 1.0)) / v0_;
     const Real rho_low = 0.0; // zero density is valid for `f` defined above.
 
-    regula_falsi(f, 0.0 /*target*/, 1.0 / _v0 /*guess*/, 
+    regula_falsi(f, 0.0 /*target*/, 1.0 / v0_ /*guess*/, 
                  rho_low /*left bracket*/, rho_high /*right bracket*/,
                  1.0e-9 /*x? tol*/, 1.0e-9 /*y? tol*/,
                  rho, &root_info, true);
@@ -310,10 +287,35 @@ class SimpleMACAW : public EosBase<SimpleMACAW> {
   static std::string EosPyType() { return EosType(); }
 
  private:
-  Real _A, _B, _Cvinf, _v0, _T0, _Gc;
+  Real A_, B_, Cvinf_, v0_, T0_, Gc_;
   MeanAtomicProperties _AZbar;
   static constexpr const unsigned long _preferred_input =
       thermalqs::density | thermalqs::specific_internal_energy;
+
+  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_);
+  }
+
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real _SieThermalPortion(
+      const Real rho, const Real T,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    const Real ratio = rho * v0_;
+    return Cvinf_ * T * (1.0 - robust::ratio(T0_, T0_ + T * std::pow(ratio, -Gc_)));
+  }
+
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real _PressureThermalPortion(
+      const Real rho, const Real T,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+      const Real numerator = rho * Cvinf_ * Gc_ * math_utils::pow<2>(T);
+      const Real denominator = T0_ * std::pow(rho * v0_, Gc_) + T;
+    return robust::ratio(numerator, denominator);
+  }
+
 };
 
 } // namespace singularity