Add lower bound for temperature for IdealFluid

Author: jyoo1042

src/PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.cpp

@@ -12,8 +12,9 @@
 
 namespace EquationsOfState {
 template <bool IsRelativistic>
-IdealFluid<IsRelativistic>::IdealFluid(const double adiabatic_index)
-    : adiabatic_index_(adiabatic_index) {}
+IdealFluid<IsRelativistic>::IdealFluid(const double adiabatic_index,
+                                       const double min_temperature)
+    : adiabatic_index_(adiabatic_index), min_temperature_(min_temperature) {}
 
 EQUATION_OF_STATE_MEMBER_DEFINITIONS(template <bool IsRelativistic>,
                                      IdealFluid<IsRelativistic>, double, 2)
@@ -44,7 +45,8 @@ bool IdealFluid<IsRelativistic>::is_equal(
 template <bool IsRelativistic>
 bool IdealFluid<IsRelativistic>::operator==(
     const IdealFluid<IsRelativistic>& rhs) const {
-  return adiabatic_index_ == rhs.adiabatic_index_;
+  return (adiabatic_index_ == rhs.adiabatic_index_) and
+         (min_temperature_ == rhs.min_temperature_);
 }
 
 template <bool IsRelativistic>
@@ -61,6 +63,7 @@ template <bool IsRelativistic>
 void IdealFluid<IsRelativistic>::pup(PUP::er& p) {
   EquationOfState<IsRelativistic, 2>::pup(p);
   p | adiabatic_index_;
+  p | min_temperature_;
 }
 
 template <bool IsRelativistic>

src/PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp

@@ -55,6 +55,14 @@ class IdealFluid : public EquationOfState<IsRelativistic, 2> {
     static constexpr Options::String help = {"Adiabatic index gamma"};
   };
 
+  struct MinTemperature {
+    using type = double;
+    static type lower_bound() { return 0.0; }
+    static constexpr Options::String help = {
+        "Minimum temperature. "
+        "This value must not be greater than or equal to 0.0"};
+  };
+
   static constexpr Options::String help = {
       "An ideal fluid equation of state.\n"
       "The pressure is related to the rest mass density by p = rho * epsilon * "
@@ -63,7 +71,7 @@ class IdealFluid : public EquationOfState<IsRelativistic, 2> {
       "index.\n"
       "The temperature T is defined as T=epsilon."};
 
-  using options = tmpl::list<AdiabaticIndex>;
+  using options = tmpl::list<AdiabaticIndex, MinTemperature>;
 
   IdealFluid() = default;
   IdealFluid(const IdealFluid&) = default;
@@ -72,7 +80,7 @@ class IdealFluid : public EquationOfState<IsRelativistic, 2> {
   IdealFluid& operator=(IdealFluid&&) = default;
   ~IdealFluid() override = default;
 
-  explicit IdealFluid(double adiabatic_index);
+  explicit IdealFluid(double adiabatic_index, double min_temperature = 0.0);
 
   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(IdealFluid, 2)
 
@@ -104,9 +112,13 @@ class IdealFluid : public EquationOfState<IsRelativistic, 2> {
 
   /// The lower bound of the specific internal energy that is valid for this EOS
   /// at the given rest mass density \f$\rho\f$
+  /// If non-zero lower bound for temperature is provided, then the lower bound
+  /// for specific internal energy is also non-zero.
   double specific_internal_energy_lower_bound(
       const double /* rest_mass_density */) const override {
-    return 0.0;
+    return (min_temperature_ == 0.0)
+               ? 0.0
+               : (min_temperature_) / (adiabatic_index_ - 1.0);
   }
 
   /// The upper bound of the specific internal energy that is valid for this EOS
@@ -119,6 +131,11 @@ class IdealFluid : public EquationOfState<IsRelativistic, 2> {
     return IsRelativistic ? 1.0 : 0.0;
   }
 
+  /// The lower bound of the temperature that is valid for this EOS.
+  /// Non-zero lower bound could be set to impose floor for the specific
+  /// internal energy.
+  double temperature_lower_bound() const override { return min_temperature_; }
+
   /// The vacuum baryon mass for this EoS
   double baryon_mass() const override {
     return hydro::units::geometric::default_baryon_mass;
@@ -128,6 +145,7 @@ class IdealFluid : public EquationOfState<IsRelativistic, 2> {
   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(2)
 
   double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
+  double min_temperature_ = std::numeric_limits<double>::signaling_NaN();
 };
 
 /// \cond

tests/Unit/PointwiseFunctions/Hydro/EquationsOfState/Test_Equilibrium3D.cpp

@@ -27,7 +27,7 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.EquationsOfState.Equilibrium3D",
   const double ideal_adiabatic_index = 2.0;
   const double cold_polytropic_index = 2.0;
   const EoS::PolytropicFluid<true> cold_eos{100.0, cold_polytropic_index};
-  const EoS::IdealFluid<true> eos_ideal_fluid{ideal_adiabatic_index};
+  const EoS::IdealFluid<true> eos_ideal_fluid{ideal_adiabatic_index, 0.0};
 
   EoS::Equilibrium3D<EoS::IdealFluid<true>> eos_ideal_fluid_3d{eos_ideal_fluid};
   const EoS::HybridEos<EoS::PolytropicFluid<true>> underlying_eos{

tests/Unit/PointwiseFunctions/Hydro/EquationsOfState/Test_IdealFluid.cpp

@@ -60,6 +60,23 @@ void check_bounds() {
         approx(hydro::units::geometric::default_baryon_mass));
 }
 
+template <bool IsRelativistic>
+void check_nonzero_bounds() {
+  const auto seed = std::random_device{}();
+  MAKE_GENERATOR(generator, seed);
+  CAPTURE(seed);
+  const auto adiabatic_index = 1.5;
+  CAPTURE(adiabatic_index);
+  auto distribution = std::uniform_real_distribution<>{1.e-9, 1.e-5};
+  auto min_temp = distribution(generator);
+  CAPTURE(min_temp);
+  const auto eos =
+      EquationsOfState::IdealFluid<IsRelativistic>{adiabatic_index, min_temp};
+  CHECK(min_temp == eos.temperature_lower_bound());
+  CHECK(min_temp / (adiabatic_index - 1) ==
+        eos.specific_internal_energy_lower_bound(1.0));
+}
+
 }  // namespace
 
 SPECTRE_TEST_CASE("Unit.PointwiseFunctions.EquationsOfState.IdealFluid",
@@ -101,29 +118,35 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.EquationsOfState.IdealFluid",
       TestHelpers::test_creation<
           std::unique_ptr<EoS::EquationOfState<true, 2>>>(
           {"IdealFluid:\n"
-           "  AdiabaticIndex: 1.6666666666666667\n"}),
+           "  AdiabaticIndex: 1.6666666666666667\n"
+           "  MinTemperature: 0.0\n"}),
       "IdealFluid", "ideal_fluid", d_for_size, 5.0 / 3.0);
   TestHelpers::EquationsOfState::check(
       TestHelpers::test_creation<
           std::unique_ptr<EoS::EquationOfState<true, 2>>>(
           {"IdealFluid:\n"
-           "  AdiabaticIndex: 1.3333333333333333\n"}),
+           "  AdiabaticIndex: 1.3333333333333333\n"
+           "  MinTemperature: 0.0\n"}),
       "IdealFluid", "ideal_fluid", dv_for_size, 4.0 / 3.0);
 
   TestHelpers::EquationsOfState::check(
       TestHelpers::test_creation<
           std::unique_ptr<EoS::EquationOfState<false, 2>>>(
           {"IdealFluid:\n"
-           "  AdiabaticIndex: 1.6666666666666667\n"}),
+           "  AdiabaticIndex: 1.6666666666666667\n"
+           "  MinTemperature: 0.0\n"}),
       "IdealFluid", "ideal_fluid", d_for_size, 5.0 / 3.0);
   TestHelpers::EquationsOfState::check(
       TestHelpers::test_creation<
           std::unique_ptr<EoS::EquationOfState<false, 2>>>(
           {"IdealFluid:\n"
-           "  AdiabaticIndex: 1.3333333333333333\n"}),
+           "  AdiabaticIndex: 1.3333333333333333\n"
+           "  MinTemperature: 0.0\n"}),
       "IdealFluid", "ideal_fluid", dv_for_size, 4.0 / 3.0);
 
   check_bounds<true>();
   check_bounds<false>();
+  check_nonzero_bounds<true>();
+  check_nonzero_bounds<false>();
   check_dominant_energy_condition_at_bound();
 }