Fix bug in Gruneisen DensityEnergyFromPressureTemperature()

CHANGELOG.md

@@ -6,6 +6,7 @@
 
 ### Fixed (Repair bugs, etc)
 - [[PR401]](https://github.com/lanl/singularity-eos/pull/401) Fix for internal energy scaling in PTE closure
+- [[PR403]](https://github.com/lanl/singularity-eos/pull/403) Fix Gruneisen EOS DensityEnergyFromPressureTemperature function
 
 ### Changed (changing behavior/API/variables/...)
 

singularity-eos/eos/eos_gruneisen.hpp

@@ -63,6 +63,8 @@ class Gruneisen : public EosBase<Gruneisen> {
         _Cv(Cv), _rho_max(RHOMAX_SAFETY * ComputeRhoMax(s1, s2, s3, rho0)) {}
   static PORTABLE_INLINE_FUNCTION Real ComputeRhoMax(const Real s1, const Real s2,
                                                      const Real s3, const Real rho0);
+  PORTABLE_INLINE_FUNCTION Real
+  MaxStableDensityAtTemperature(const Real temperature) const;
   Gruneisen GetOnDevice() { return *this; }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy(
@@ -397,28 +399,87 @@ PORTABLE_INLINE_FUNCTION Real Gruneisen::BulkModulusFromDensityTemperature(
   return BulkModulusFromDensityInternalEnergy(
       rho, InternalEnergyFromDensityTemperature(rho, temp));
 }
+PORTABLE_INLINE_FUNCTION Real
+Gruneisen::MaxStableDensityAtTemperature(const Real temperature) const {
+  // Because of the constant cold curve assumption, this EOS is thermodynamically
+  // inconsistent, and leads to states that are effectively off the EOS surface even
+  // though the temperature is positive. Mathematically, this means that the \Gamma\rho(e
+  // - e_H) term becomes highly negative and dominates the positive P_H term at some
+  // density. This results in a local maximum in the pressure as a function of density for
+  // a given temperature. Beyond this point, the pressure decreases with increasing
+  // density, which is thermodynamcially unstable. In a thermodynamically consistent EOS,
+  // the cold curve energy would increase with density, leading an appropriate bounding at
+  // T=0.
+
+  // Since E_H and P_H are monotonic up to the singularity given by _rho_max, if the
+  // derivative of the pressure is negative at _rho_max, a maximum exists and this should
+  // be the highest density for the isotherm. If this maximum doesn't exist, then the EOS
+  // is thermodynamically consistent up to the maximum density for this isotherm.
+  Real slope_at_max_density =
+      dPres_drho_e(_rho_max, InternalEnergyFromDensityTemperature(_rho_max, temperature));
+  if (slope_at_max_density >= 0) {
+    // No maximum pressure before _rho_max
+    return _rho_max;
+  }
+
+  // Maximum pressure should exist... do a root find to locate where the derivative is
+  // zero
+  auto dPdrho_T = PORTABLE_LAMBDA(const Real r) {
+    return dPres_drho_e(r, InternalEnergyFromDensityTemperature(r, temperature));
+  };
+  Real rho_lower = 0.9 * _rho0;
+  Real rho_upper = std::min(_rho_max, 1.0e4);
+  Real rho_at_max_P;
+  using RootFinding1D::regula_falsi;
+  using RootFinding1D::Status;
+  auto status = regula_falsi(dPdrho_T, 0., _rho0, rho_lower, rho_upper, 1.0e-8, 1.0e-8,
+                             rho_at_max_P);
+  if (status != Status::SUCCESS) {
+    // Root finder failed even though the solution should be bracketed
+    EOS_ERROR("Gruneisen::MaxStableDensityAtTemperature: "
+              "Root find failed to find maximum P at T despite aparent bracket\n");
+  }
+  return rho_at_max_P;
+}
 template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature(
     const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const {
-  sie = _Cv * (temp - _T0);
   // We have a branch at rho0, so we need to decide, based on our pressure, whether we
   // should be above or below rho0
-  Real Pref = PressureFromDensityInternalEnergy(_rho0, sie);
+  Real Pref = PressureFromDensityTemperature(_rho0, temp);
+  Real rho_lower;
+  Real rho_upper;
+  // Pick bounds appropriate depending on whether in compression or expansion
   if (press < Pref) {
-    rho = (press - _P0 + _C0 * _C0 * _rho0) / (_C0 * _C0 + _G0 * sie);
-  } else { // We are in compression; iterate
-    auto PofRatE = [&](const Real r) {
-      return PressureFromDensityInternalEnergy(r, sie);
-    };
-    using RootFinding1D::regula_falsi;
-    using RootFinding1D::Status;
-    auto status = regula_falsi(PofRatE, press, _rho0, 1.0e-5, 1.0e3, 1.0e-8, 1.0e-8, rho);
-    if (status != Status::SUCCESS) {
-      // Root finder failed even though the solution was bracketed... this is an error
-      EOS_ERROR("Gruneisen::DensityEnergyFromPressureTemperature: "
-                "Root find failed to find a solution given P, T\n");
+    rho_lower = 0.;
+    rho_upper = 1.1 * _rho0;
+  } else {
+    rho_lower = 0.9 * _rho0;
+    // Find maximum thermodynamically _stable_ density at this temperature. Use this to
+    // check if we're actually on the EOS surface or not
+    rho_upper = MaxStableDensityAtTemperature(temp);
+    auto pres_max = PressureFromDensityTemperature(rho_upper, temp);
+    if (press > pres_max) {
+      // We're off the EOS surface
+      using PortsOfCall::printf;
+      printf("Requested pressure, %.15g, exceeds maximum, %.15g, for temperature, %.15g",
+             press, pres_max, temp);
+      PORTABLE_ALWAYS_THROW_OR_ABORT("Input pressure is off EOS surface");
     }
   }
+  auto PofRatT = PORTABLE_LAMBDA(const Real r) {
+    return PressureFromDensityTemperature(r, temp);
+  };
+  using RootFinding1D::regula_falsi;
+  using RootFinding1D::Status;
+  auto status =
+      regula_falsi(PofRatT, press, _rho0, rho_lower, rho_upper, 1.0e-8, 1.0e-8, rho);
+  if (status != Status::SUCCESS) {
+    // Root finder failed even though the solution was bracketed... this is an error
+    EOS_ERROR("Gruneisen::DensityEnergyFromPressureTemperature: "
+              "Root find failed to find a solution given P, T\n");
+  }
+  sie = InternalEnergyFromDensityTemperature(rho, temp);
 }
 template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION void

test/test_eos_gruneisen.cpp

@@ -1,5 +1,5 @@
 //------------------------------------------------------------------------------
-// © 2021-2023. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract 89233218CNA000001
 // for Los Alamos National Laboratory (LANL), which is operated by Triad
 // National Security, LLC for the U.S.  Department of Energy/National
@@ -33,6 +33,8 @@ PORTABLE_INLINE_FUNCTION Real QuadFormulaMinus(Real a, Real b, Real c) {
   return (-b - std::sqrt(b * b - 4 * a * c)) / (2 * a);
 }
 
+constexpr Real REAL_TOL = std::numeric_limits<Real>::epsilon() * 1e4; // 1e-12 for double
+
 // Only run exception checking when we aren't offloading to the GPUs
 #ifdef PORTABILITY_STRATEGY_NONE
 SCENARIO("Gruneisen EOS entropy is disabled", "[GruneisenEOS][Entropy]") {
@@ -206,40 +208,57 @@ SCENARIO("Gruneisen EOS", "[VectorEOS][GruneisenEOS]") {
   }
 }
 
-SCENARIO("Aluminum Gruneisen EOS Sound Speed and Pressure Comparison", "[GruneisenEOS]") {
-  GIVEN("Parameters for a Gruneisen EOS") {
+SCENARIO("Aluminum Gruneisen EOS", "[GruneisenEOS]") {
+  GIVEN("Parameters for an aluminum Gruneisen EOS") {
     // Unit conversions
     // constexpr Real mm = 10.;
     constexpr Real cm = 1.;
     constexpr Real us = 1.e-06;
     constexpr Real Mbar = 1.e12;
     constexpr Real Mbcc_per_g = 1e12;
-    // Gruneisen parameters for copper
-    constexpr Real C0 = 0.535 * cm / us;
-    constexpr Real S1 = 1.34;
+    // Gruneisen parameters for aluminum
+    constexpr Real C0 = 0.524 * cm / us;
+    constexpr Real S1 = 1.4;
     constexpr Real S2 = 0.;
     constexpr Real S3 = 0.;
     constexpr Real Gamma0 = 1.97;
-    constexpr Real b = 0.;
-    constexpr Real rho0 = 2.714000;
+    constexpr Real b = 0.48;
+    constexpr Real rho0 = 2.703;
     constexpr Real T0 = 298.;
-    constexpr Real P0 = 1e-06 * Mbar;
+    constexpr Real P0 = 0. * Mbar;
     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();
+    GIVEN("Ambient pressure and temperature") {
+      constexpr Real P_ambient = 1.0e-06 * Mbar;
+      constexpr Real T_ambient = 298.;
+      WHEN("A DensityEnergyFromPressureTemperature() lookup is performed") {
+        Real test_density;
+        Real test_energy;
+        Real *lambda;
+        THEN("The root find should converge") {
+          eos.DensityEnergyFromPressureTemperature(P_ambient, T_ambient, lambda,
+                                                   test_density, test_energy);
+        }
+      }
+    }
     GIVEN("Density and energy") {
-      constexpr Real density = 5.92418956756592;            // g/cm^3
-      constexpr Real energy = 792486007.804619;             // erg/g
-      constexpr Real true_pres = 2.620656373250729;         // Mbar
-      constexpr Real true_sound_speed = 1.5247992468363685; // cm/us
+      constexpr Real density = 5.92418956756592;          // g/cm^3
+      constexpr Real energy = 792486007.804619;           // erg/g
+      constexpr Real true_pres = 2.19200396047868;        // Mbar
+      constexpr Real true_sound_speed = 1.29592658233752; // cm/us
+      THEN("The density should not exceeed the computed rho_max") {
+        const Real density_max = Gruneisen::ComputeRhoMax(S1, S2, S3, rho0);
+        REQUIRE(density < density_max);
+      }
       WHEN("A P(rho, e) lookup is performed") {
         Real pres = eos.PressureFromDensityInternalEnergy(density, energy);
         THEN("The correct pressure should be returned") {
           pres = pres / Mbar;
           INFO("Density: " << density << "  Energy: " << energy << "  Pressure: " << pres
                            << " Mbar  True pressure: " << true_pres << " Mbar");
-          REQUIRE(isClose(pres, true_pres, 1e-12));
+          REQUIRE(isClose(pres, true_pres, REAL_TOL));
         }
       }
       WHEN("A B_S(rho, e) lookup is performed") {
@@ -251,7 +270,29 @@ SCENARIO("Aluminum Gruneisen EOS Sound Speed and Pressure Comparison", "[Gruneis
                            << "  Sound speed: " << sound_speed
                            << " cm/us  True sound speed: " << true_sound_speed
                            << " cm/us");
-          REQUIRE(isClose(sound_speed, true_sound_speed, 1e-12));
+          REQUIRE(isClose(sound_speed, true_sound_speed, REAL_TOL));
+        }
+      }
+      WHEN("The pressure and temperature are determined from the density and energy") {
+        const Real temperature =
+            eos.TemperatureFromDensityInternalEnergy(density, energy);
+        const Real pressure = eos.PressureFromDensityInternalEnergy(density, energy);
+        AND_WHEN("A DensityEnergyFromPressureTemperature() lookup is performed") {
+          Real test_density;
+          Real test_energy;
+          Real *lambda;
+          eos.DensityEnergyFromPressureTemperature(pressure, temperature, lambda,
+                                                   test_density, test_energy);
+          THEN("The consistent energy and density should be returned") {
+            INFO("Pressure:     " << pressure << " microbar"
+                                  << "  Temperature: " << temperature << " K       ");
+            INFO("Density:      " << density << " g/cm^3      "
+                                  << "  Energy:      " << energy << " erg/g   ");
+            INFO("Calc Density: " << test_density << " g/cm^3      "
+                                  << "  Calc Energy: " << test_energy << " erg/g   ");
+            CHECK(isClose(density, test_density, REAL_TOL));
+            CHECK(isClose(energy, test_energy, REAL_TOL));
+          }
         }
       }
       WHEN("A finite difference approximation is used for the bulk modulus") {
@@ -294,7 +335,7 @@ SCENARIO("Aluminum Gruneisen EOS Sound Speed and Pressure Comparison", "[Gruneis
             INFO("Density: " << density << "  Energy: " << energy
                              << "  Pressure: " << pres / Mbar
                              << " Mbar  True pressure: " << true_pres / Mbar << " Mbar");
-            REQUIRE(isClose(pres, true_pres, 1e-12));
+            REQUIRE(isClose(pres, true_pres, REAL_TOL));
           }
         }
       }
@@ -399,7 +440,7 @@ SCENARIO("Gruneisen EOS density limit") {
         THEN("The generated rho_max parameter should be properly set") {
           const Real rho_max = eos.ComputeRhoMax(S1, S2, S3, rho0);
           INFO("True rho_max: " << rho_max_true << ", Calculated rho_max:" << rho_max);
-          REQUIRE(isClose(rho_max, rho_max_true, 1e-12));
+          REQUIRE(isClose(rho_max, rho_max_true, REAL_TOL));
         }
         WHEN("Lookups are performed beyond the maximum density") {
           // Note: there is a small safety factor to prevent us from hitting the
@@ -495,7 +536,7 @@ SCENARIO("Gruneisen EOS density limit") {
             const Real rho_max_true = rho0 / (1 - eta_max);
             const Real rho_max = eos.ComputeRhoMax(S1, S2, S3, rho0);
             INFO("True rho_max: " << rho_max_true << ", Calculated rho_max:" << rho_max);
-            REQUIRE(isClose(rho_max, rho_max_true, 1e-12));
+            REQUIRE(isClose(rho_max, rho_max_true, REAL_TOL));
           }
         }
         WHEN("The quadratic oot multiplicity is 2") {
@@ -511,7 +552,7 @@ SCENARIO("Gruneisen EOS density limit") {
             const Real rho_max_true = rho0 / (1 - eta_max);
             const Real rho_max = eos.ComputeRhoMax(S1, S2, S3, rho0);
             INFO("True rho_max: " << rho_max_true << ", Calculated rho_max:" << rho_max);
-            REQUIRE(isClose(rho_max, rho_max_true, 1e-12));
+            REQUIRE(isClose(rho_max, rho_max_true, REAL_TOL));
           }
         }
         WHEN("No root exists") {