Add protection against negative pressure in HybridEos

src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.tpp

@@ -140,9 +140,11 @@ class FunctionOfMu {
       eps_min = equation_of_state_.specific_internal_energy_lower_bound(
           rest_mass_density_times_lorentz_factor_ / lorentz_max,
           electron_fraction_);
-    } else {
+    } else if constexpr (EosType::thermodynamic_dim == 2) {
       eps_min = equation_of_state_.specific_internal_energy_lower_bound(
           rest_mass_density_times_lorentz_factor_ / lorentz_max);
+    } else {
+      eps_min = equation_of_state_.specific_internal_energy_lower_bound();
     }
     q_ = tau / rest_mass_density_times_lorentz_factor;
     if constexpr (EnforcePhysicality) {
@@ -295,11 +297,16 @@ Primitives FunctionOfMu<EnforcePhysicality, EosType>::primitives(
                        rho_hat, electron_fraction_),
                    equation_of_state_.specific_internal_energy_upper_bound(
                        rho_hat, electron_fraction_));
-  } else {
+  } else if constexpr (EosType::thermodynamic_dim == 2) {
     epsilon_hat = std::clamp(
         epsilon_hat,
         equation_of_state_.specific_internal_energy_lower_bound(rho_hat),
         equation_of_state_.specific_internal_energy_upper_bound(rho_hat));
+  } else {
+    epsilon_hat = std::clamp(
+        epsilon_hat,
+        equation_of_state_.specific_internal_energy_lower_bound(),
+        equation_of_state_.specific_internal_energy_upper_bound());
   }
   // Pressure from EOS
   double p_hat = std::numeric_limits<double>::signaling_NaN();

src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAlHydro.tpp

@@ -57,9 +57,12 @@ class FunctionOfZ {
     if constexpr (EosType::thermodynamic_dim == 3) {
       eps_min = equation_of_state_.specific_internal_energy_lower_bound(
           rho_min, electron_fraction);
-    } else {
+    } else if constexpr (EosType::thermodynamic_dim == 2) {
       eps_min =
           equation_of_state_.specific_internal_energy_lower_bound(rho_min);
+    } else {
+      eps_min =
+          equation_of_state_.specific_internal_energy_lower_bound();
     }
 
     if constexpr (EnforcePhysicality) {
@@ -134,11 +137,16 @@ Primitives FunctionOfZ<EosType, EnforcePhysicality>::primitives(
                        rho_hat, electron_fraction_),
                    equation_of_state_.specific_internal_energy_upper_bound(
                        rho_hat, electron_fraction_));
-  } else {
+  } else if constexpr (EosType::thermodynamic_dim == 2) {
     epsilon_hat = std::clamp(
         epsilon_hat,
         equation_of_state_.specific_internal_energy_lower_bound(rho_hat),
         equation_of_state_.specific_internal_energy_upper_bound(rho_hat));
+  } else {
+    epsilon_hat = std::clamp(
+        epsilon_hat,
+        equation_of_state_.specific_internal_energy_lower_bound(),
+        equation_of_state_.specific_internal_energy_upper_bound());
   }
 
   // Pressure from EOS

src/PointwiseFunctions/Hydro/EquationsOfState/Barotropic2D.hpp

@@ -142,17 +142,15 @@ class Barotropic2D : public EquationOfState<ColdEos::is_relativistic, 2> {
   /// The lower bound of the specific internal energy that is valid for this EOS
   /// at the given rest mass density \f$\rho\f$ and electron fraction \f$Y_e\f$
   double specific_internal_energy_lower_bound(
-      const double rest_mass_density) const override {
-    return underlying_eos_.specific_internal_energy_lower_bound(
-        rest_mass_density);
+      const double /*rest_mass_density*/) const override {
+    return underlying_eos_.specific_internal_energy_lower_bound();
   }
 
   /// The upper bound of the specific internal energy that is valid for this EOS
   /// at the given rest mass density \f$\rho\f$
   double specific_internal_energy_upper_bound(
-      const double rest_mass_density) const override {
-    return underlying_eos_.specific_internal_energy_upper_bound(
-        rest_mass_density);
+      const double /*rest_mass_density*/) const override {
+    return underlying_eos_.specific_internal_energy_upper_bound();
   }
 
   /// The lower bound of the specific enthalpy that is valid for this EOS

src/PointwiseFunctions/Hydro/EquationsOfState/Barotropic3D.hpp

@@ -147,19 +147,17 @@ class Barotropic3D : public EquationOfState<ColdEquilEos::is_relativistic, 3> {
   /// The lower bound of the specific internal energy that is valid for this EOS
   /// at the given rest mass density \f$\rho\f$ and electron fraction \f$Y_e\f$
   double specific_internal_energy_lower_bound(
-      const double rest_mass_density,
+      const double /*rest_mass_density*/,
       const double /*electron_fraction*/) const override {
-    return underlying_eos_.specific_internal_energy_lower_bound(
-        rest_mass_density);
+    return underlying_eos_.specific_internal_energy_lower_bound();
   }
 
   /// The upper bound of the specific internal energy that is valid for this EOS
   /// at the given rest mass density \f$\rho\f$
   double specific_internal_energy_upper_bound(
-      const double rest_mass_density,
+      const double /*rest_mass_density*/,
       const double /*electron_fraction*/) const override {
-    return underlying_eos_.specific_internal_energy_upper_bound(
-        rest_mass_density);
+    return underlying_eos_.specific_internal_energy_upper_bound();
   }
 
   /// The lower bound of the specific enthalpy that is valid for this EOS

src/PointwiseFunctions/Hydro/EquationsOfState/Enthalpy.hpp

@@ -248,23 +248,17 @@ class Enthalpy : public EquationOfState<true, 1> {
     return std::numeric_limits<double>::max();
   }
 
+  /// The lower bound of the specific enthalpy that is valid for this EOS
+  double specific_enthalpy_lower_bound() const override { return 1.0; }
+
   /// The lower bound of the specific internal energy that is valid for this EOS
-  /// at the given rest mass density \f$\rho\f$
-  double specific_internal_energy_lower_bound(
-      const double /* rest_mass_density */) const override {
-    return 0.0;
-  }
+  double specific_internal_energy_lower_bound() const override { return 0.0; }
 
   /// The upper bound of the specific internal energy that is valid for this EOS
-  /// at the given rest mass density \f$\rho\f$
-  double specific_internal_energy_upper_bound(
-      const double /* rest_mass_density */) const override {
+  double specific_internal_energy_upper_bound() const override {
     return std::numeric_limits<double>::max();
   }
 
-  /// The lower bound of the specific enthalpy that is valid for this EOS
-  double specific_enthalpy_lower_bound() const override { return 1.0; }
-
   /// The vacuum baryon mass for this EoS
   double baryon_mass() const override { return low_density_eos_.baryon_mass(); }
 

src/PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp

@@ -308,15 +308,14 @@ class EquationOfState<IsRelativistic, 1> : public PUP::able {
   virtual double temperature_upper_bound() const {
     return std::numeric_limits<double>::max();
   };
+
   /// The lower bound of the specific internal energy that is valid for this EOS
-  /// at the given rest mass density \f$\rho\f$
-  virtual double specific_internal_energy_lower_bound(
-      double rest_mass_density) const = 0;
+  virtual double specific_internal_energy_lower_bound() const { return 0.0; };
 
   /// The upper bound of the specific internal energy that is valid for this EOS
-  /// at the given rest mass density \f$\rho\f$
-  virtual double specific_internal_energy_upper_bound(
-      double rest_mass_density) const = 0;
+  virtual double specific_internal_energy_upper_bound() const {
+    return std::numeric_limits<double>::max();
+  };
 
   /// The lower bound of the specific enthalpy that is valid for this EOS
   virtual double specific_enthalpy_lower_bound() const = 0;

src/PointwiseFunctions/Hydro/EquationsOfState/HybridEos.cpp

@@ -81,12 +81,14 @@ Scalar<DataType>
 HybridEos<ColdEquationOfState>::pressure_from_density_and_energy_impl(
     const Scalar<DataType>& rest_mass_density,
     const Scalar<DataType>& specific_internal_energy) const {
+  using std::max;
   return Scalar<DataType>{
       get(cold_eos_.pressure_from_density(rest_mass_density)) +
       get(rest_mass_density) * (thermal_adiabatic_index_ - 1.0) *
-          (get(specific_internal_energy) -
-           get(cold_eos_.specific_internal_energy_from_density(
-               rest_mass_density)))};
+          max((get(specific_internal_energy) -
+               get(cold_eos_.specific_internal_energy_from_density(
+                   rest_mass_density))),
+              0.0)};
 }
 
 template <typename ColdEquationOfState>
@@ -95,21 +97,22 @@ Scalar<DataType>
 HybridEos<ColdEquationOfState>::pressure_from_density_and_enthalpy_impl(
     const Scalar<DataType>& rest_mass_density,
     const Scalar<DataType>& specific_enthalpy) const {
+  using std::max;
   if constexpr (ColdEquationOfState::is_relativistic) {
     return Scalar<DataType>{
         (get(cold_eos_.pressure_from_density(rest_mass_density)) +
          get(rest_mass_density) * (thermal_adiabatic_index_ - 1.0) *
-             (get(specific_enthalpy) - 1.0 -
-              get(cold_eos_.specific_internal_energy_from_density(
-                  rest_mass_density)))) /
+             max((get(specific_enthalpy) - 1.0 -
+                  get(cold_eos_.specific_internal_energy_from_density(
+                      rest_mass_density))), 0.0)) /
         thermal_adiabatic_index_};
   } else {
     return Scalar<DataType>{
         (get(cold_eos_.pressure_from_density(rest_mass_density)) +
          get(rest_mass_density) * (thermal_adiabatic_index_ - 1.0) *
-             (get(specific_enthalpy) -
-              get(cold_eos_.specific_internal_energy_from_density(
-                  rest_mass_density)))) /
+             max((get(specific_enthalpy) -
+                  get(cold_eos_.specific_internal_energy_from_density(
+                      rest_mass_density))), 0.0)) /
         thermal_adiabatic_index_};
   }
 }
@@ -120,11 +123,12 @@ Scalar<DataType> HybridEos<ColdEquationOfState>::
     specific_internal_energy_from_density_and_pressure_impl(
         const Scalar<DataType>& rest_mass_density,
         const Scalar<DataType>& pressure) const {
+  using std::max;
   return Scalar<DataType>{
       get(cold_eos_.specific_internal_energy_from_density(rest_mass_density)) +
       1.0 / (thermal_adiabatic_index_ - 1.0) *
-          (get(pressure) -
-           get(cold_eos_.pressure_from_density(rest_mass_density))) /
+          max((get(pressure) -
+               get(cold_eos_.pressure_from_density(rest_mass_density))), 0.0) /
           get(rest_mass_density)};
 }
 
@@ -134,10 +138,11 @@ Scalar<DataType>
 HybridEos<ColdEquationOfState>::temperature_from_density_and_energy_impl(
     const Scalar<DataType>& rest_mass_density,
     const Scalar<DataType>& specific_internal_energy) const {
+  using std::max;
   return Scalar<DataType>{(thermal_adiabatic_index_ - 1.0) *
-                          (get(specific_internal_energy) -
+                          max((get(specific_internal_energy) -
                            get(cold_eos_.specific_internal_energy_from_density(
-                               rest_mass_density)))};
+                               rest_mass_density))), 0.0)};
 }
 
 template <typename ColdEquationOfState>
@@ -173,14 +178,15 @@ Scalar<DataType> HybridEos<ColdEquationOfState>::
     kappa_times_p_over_rho_squared_from_density_and_energy_impl(
         const Scalar<DataType>& rest_mass_density,
         const Scalar<DataType>& specific_internal_energy) const {
+  using std::max;
   return Scalar<DataType>{
       (thermal_adiabatic_index_ - 1.0) *
           get(cold_eos_.pressure_from_density(rest_mass_density)) /
           get(rest_mass_density) +
       square(thermal_adiabatic_index_ - 1.0) *
-          (get(specific_internal_energy) -
-           get(cold_eos_.specific_internal_energy_from_density(
-               rest_mass_density)))};
+          max((get(specific_internal_energy) -
+               get(cold_eos_.specific_internal_energy_from_density(
+                   rest_mass_density))), 0.0)};
 }
 }  // namespace EquationsOfState
 

src/PointwiseFunctions/Hydro/EquationsOfState/HybridEos.hpp

@@ -141,7 +141,8 @@ class HybridEos
   /// at the given rest mass density \f$\rho\f$
   double specific_internal_energy_lower_bound(
       const double rest_mass_density) const override {
-    return cold_eos_.specific_internal_energy_lower_bound(rest_mass_density);
+    return get(cold_eos_.specific_internal_energy_from_density(
+        Scalar<double>{rest_mass_density}));
   }
 
   /// The upper bound of the specific internal energy that is valid for this EOS

src/PointwiseFunctions/Hydro/EquationsOfState/PiecewisePolytropicFluid.cpp

@@ -302,6 +302,24 @@ double PiecewisePolytropicFluid<IsRelativistic>::rest_mass_density_upper_bound()
   }
   return std::numeric_limits<double>::max();
 }
+
+template <bool IsRelativistic>
+double PiecewisePolytropicFluid<
+    IsRelativistic>::specific_internal_energy_upper_bound() const {
+  // this bound comes from the dominant energy condition which implies
+  // that the pressure is bounded by the total energy density,
+  // i.e. p < e = rho * (1 + eps)
+  if (IsRelativistic and polytropic_exponent_hi_ > 2.0) {
+    const double eint_boundary_constant =
+        (polytropic_exponent_hi_ - polytropic_exponent_lo_) *
+        polytropic_constant_lo_ /
+        ((polytropic_exponent_hi_ - 1.0) * (polytropic_exponent_lo_ - 1.0)) *
+        pow(transition_density_, polytropic_exponent_lo_ - 1.0);
+    return (1.0 + (polytropic_exponent_hi_ - 1.0) * eint_boundary_constant) /
+           (polytropic_exponent_hi_ - 2.0);
+  }
+  return std::numeric_limits<double>::max();
+}
 }  // namespace EquationsOfState
 
 template class EquationsOfState::PiecewisePolytropicFluid<true>;

src/PointwiseFunctions/Hydro/EquationsOfState/PiecewisePolytropicFluid.hpp

@@ -145,25 +145,17 @@ class PiecewisePolytropicFluid : public EquationOfState<IsRelativistic, 1> {
   /// The upper bound of the rest mass density that is valid for this EOS
   double rest_mass_density_upper_bound() const override;
 
-  /// The lower bound of the specific internal energy that is valid for this EOS
-  /// at the given rest mass density \f$\rho\f$
-  double specific_internal_energy_lower_bound(
-      const double /* rest_mass_density */) const override {
-    return 0.0;
-  }
-
-  /// The upper bound of the specific internal energy that is valid for this EOS
-  /// at the given rest mass density \f$\rho\f$
-  double specific_internal_energy_upper_bound(
-      const double /* rest_mass_density */) const override {
-    return std::numeric_limits<double>::max();
-  }
-
   /// The lower bound of the specific enthalpy that is valid for this EOS
   double specific_enthalpy_lower_bound() const override {
     return IsRelativistic ? 1.0 : 0.0;
   }
 
+  /// The lower bound of the specific internal energy that is valid for this EOS
+  double specific_internal_energy_lower_bound() const override { return 0.0; }
+
+  /// The upper bound of the specific internal energy that is valid for this EOS
+  double specific_internal_energy_upper_bound() const override;
+
   /// The vacuum baryon mass for this EoS
   double baryon_mass() const override {
     return hydro::units::geometric::default_baryon_mass;

src/PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.cpp

@@ -144,6 +144,18 @@ double PolytropicFluid<IsRelativistic>::rest_mass_density_upper_bound() const {
   }
   return std::numeric_limits<double>::max();
 }
+
+template <bool IsRelativistic>
+double PolytropicFluid<IsRelativistic>::specific_internal_energy_upper_bound()
+    const {
+  // this bound comes from the dominant energy condition which implies
+  // that the pressure is bounded by the total energy density,
+  // i.e. p < e = rho * (1 + eps)
+  if (IsRelativistic and polytropic_exponent_ > 2.0) {
+    return 1.0 / (polytropic_exponent_ - 2.0);
+  }
+  return std::numeric_limits<double>::max();
+}
 }  // namespace EquationsOfState
 
 template class EquationsOfState::PolytropicFluid<true>;

src/PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp

@@ -104,25 +104,17 @@ class PolytropicFluid : public EquationOfState<IsRelativistic, 1> {
   /// The upper bound of the rest mass density that is valid for this EOS
   double rest_mass_density_upper_bound() const override;
 
-  /// The lower bound of the specific internal energy that is valid for this EOS
-  /// at the given rest mass density \f$\rho\f$
-  double specific_internal_energy_lower_bound(
-      const double /* rest_mass_density */) const override {
-    return 0.0;
-  }
-
-  /// The upper bound of the specific internal energy that is valid for this EOS
-  /// at the given rest mass density \f$\rho\f$
-  double specific_internal_energy_upper_bound(
-      const double /* rest_mass_density */) const override {
-    return std::numeric_limits<double>::max();
-  }
-
   /// The lower bound of the specific enthalpy that is valid for this EOS
   double specific_enthalpy_lower_bound() const override {
     return IsRelativistic ? 1.0 : 0.0;
   }
 
+  /// The lower bound of the specific internal energy that is valid for this EOS
+  double specific_internal_energy_lower_bound() const override { return 0.0; }
+
+  /// The upper bound of the specific internal energy that is valid for this EOS
+  double specific_internal_energy_upper_bound() const override;
+
   /// The vacuum baryon mass for this EoS
   double baryon_mass() const override {
     return hydro::units::geometric::default_baryon_mass;