Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <boost/preprocessor/arithmetic/dec.hpp>
       7             : #include <boost/preprocessor/arithmetic/inc.hpp>
       8             : #include <boost/preprocessor/control/expr_iif.hpp>
       9             : #include <boost/preprocessor/list/adt.hpp>
      10             : #include <boost/preprocessor/repetition/for.hpp>
      11             : #include <boost/preprocessor/repetition/repeat.hpp>
      12             : #include <boost/preprocessor/tuple/to_list.hpp>
      13             : #include <cstddef>
      14             : #include <limits>
      15             : #include <pup.h>
      16             : 
      17             : #include "DataStructures/Tensor/TypeAliases.hpp"
      18             : #include "Options/String.hpp"
      19             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      20             : #include "PointwiseFunctions/Hydro/Units.hpp"
      21             : #include "Utilities/Serialization/CharmPupable.hpp"
      22             : #include "Utilities/TMPL.hpp"
      23             : 
      24             : /// \cond
      25             : class DataVector;
      26             : /// \endcond
      27             : 
      28             : namespace EquationsOfState {
      29             : /*!
      30             :  * \ingroup EquationsOfStateGroup
      31             :  *
      32             :  * \brief Hybrid equation of state combining a barotropic EOS for cold
      33             :  * (zero-temperature) part with a simple thermal part
      34             :  *
      35             :  * The hybrid equation of state:
      36             :  *
      37             :  * \f[
      38             :  * p = p_{cold}(\rho) + \rho (\Gamma_{th}-1) (\epsilon - \epsilon_{cold}(\rho))
      39             :  * \f]
      40             :  *
      41             :  * where \f$p\f$ is the pressure, \f$\rho\f$ is the rest mass density,
      42             :  * \f$\epsilon\f$ is the specific internal energy, \f$p_{cold}\f$ and
      43             :  * \f$\epsilon_{cold}\f$ are the pressure and specific internal energy evaluated
      44             :  * using the cold EOS, and \f$\Gamma_{th}\f$ is the adiabatic index for the
      45             :  * thermal part.
      46             :  *
      47             :  * The temperature \f$T\f$ is defined as
      48             :  *
      49             :  * \f[
      50             :  * T = (\Gamma_{th} - 1) (\epsilon - \epsilon_{cold})
      51             :  * \f]
      52             :  *
      53             :  * This is the amount of internal energy above that of the cold EOS.
      54             :  *
      55             :  * For a hybrid EOS with the cold EOS being a polytrope we have
      56             :  *
      57             :  * \f{align}
      58             :  * p &= \kappa \rho^\Gamma + \rho T
      59             :  *     =\kappa\rho^\Gamma+\rho
      60             :  *     (\Gamma-1)\left(\epsilon-\frac{K\rho^{\Gamma-1}}{\Gamma-1}\right), \\
      61             :  * T&=(\Gamma-1)\left(\epsilon-\frac{K\rho^{\Gamma-1}}{\Gamma-1}\right), \\
      62             :  * \epsilon &= \frac{1}{(\Gamma-1)}\frac{p}{\rho}
      63             :  * \f}
      64             :  */
      65             : template <typename ColdEquationOfState>
      66           1 : class HybridEos
      67             :     : public EquationOfState<ColdEquationOfState::is_relativistic, 2> {
      68             :  public:
      69           0 :   static constexpr size_t thermodynamic_dim = 2;
      70           0 :   static constexpr bool is_relativistic = ColdEquationOfState::is_relativistic;
      71             : 
      72           0 :   struct ColdEos {
      73           0 :     using type = ColdEquationOfState;
      74           0 :     static constexpr Options::String help = {"Cold equation of state"};
      75           0 :     static std::string name() {
      76             :       return pretty_type::short_name<ColdEquationOfState>();
      77             :     }
      78             :   };
      79             : 
      80           0 :   struct ThermalAdiabaticIndex {
      81           0 :     using type = double;
      82           0 :     static constexpr Options::String help = {"Adiabatic index Gamma_th"};
      83             :   };
      84             : 
      85           0 :   struct MinTemperature {
      86           0 :     using type = double;
      87           0 :     static type lower_bound() { return 0.0; }
      88           0 :     static constexpr Options::String help = {
      89             :         "Minimum temperature. "
      90             :         "This value must be non-negative."};
      91             :   };
      92             : 
      93           0 :   static constexpr Options::String help = {
      94             :       "A hybrid equation of state combining a cold EOS with a simple thermal "
      95             :       "part.  The pressure is related to the rest mass density by "
      96             :       " p = p_cold(rho) + rho * (Gamma_th - 1) * (epsilon - "
      97             :       "epsilon_cold(rho)), where p is the pressure, rho is the rest mass "
      98             :       "density, epsilon is the specific internal energy, p_cold and "
      99             :       "epsilon_cold are the pressure and specific internal energy evaluated "
     100             :       "using the cold EOS and Gamma_th is the adiabatic index for the thermal "
     101             :       "part."};
     102             : 
     103           0 :   using options = tmpl::list<ColdEos, ThermalAdiabaticIndex, MinTemperature>;
     104             : 
     105           0 :   HybridEos() = default;
     106           0 :   HybridEos(const HybridEos&) = default;
     107           0 :   HybridEos& operator=(const HybridEos&) = default;
     108           0 :   HybridEos(HybridEos&&) = default;
     109           0 :   HybridEos& operator=(HybridEos&&) = default;
     110           0 :   ~HybridEos() override = default;
     111             : 
     112           0 :   HybridEos(ColdEquationOfState cold_eos, double thermal_adiabatic_index,
     113             :             double min_temperature = 0.0);
     114             : 
     115             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(HybridEos, 2)
     116             : 
     117           0 :   WRAPPED_PUPable_decl_base_template(  // NOLINT
     118             :       SINGLE_ARG(EquationOfState<is_relativistic, 2>), HybridEos);
     119             : 
     120           0 :   std::unique_ptr<EquationOfState<is_relativistic, 2>> get_clone()
     121             :       const override;
     122             : 
     123           0 :   std::unique_ptr<EquationOfState<is_relativistic, 3>> promote_to_3d_eos()
     124             :       const override;
     125             : 
     126             :   /// \brief Returns `true` if the EOS is barotropic
     127           1 :   bool is_barotropic() const override { return false; }
     128             : 
     129           0 :   bool operator==(const HybridEos<ColdEquationOfState>& rhs) const;
     130             : 
     131           0 :   bool operator!=(const HybridEos<ColdEquationOfState>& rhs) const;
     132             : 
     133           0 :   bool is_equal(const EquationOfState<is_relativistic, 2>& rhs) const override;
     134             : 
     135           0 :   static std::string name() {
     136             :     return "HybridEos(" + pretty_type::name<ColdEquationOfState>() + ")";
     137             :   }
     138             : 
     139             :   /// The lower bound of the rest mass density that is valid for this EOS
     140           1 :   double rest_mass_density_lower_bound() const override {
     141             :     return cold_eos_.rest_mass_density_lower_bound();
     142             :   }
     143             : 
     144             :   /// The upper bound of the rest mass density that is valid for this EOS
     145           1 :   double rest_mass_density_upper_bound() const override {
     146             :     return cold_eos_.rest_mass_density_upper_bound();
     147             :   }
     148             : 
     149             :   /// The lower bound of the specific internal energy that is valid for this EOS
     150             :   /// at the given rest mass density \f$\rho\f$
     151           1 :   double specific_internal_energy_lower_bound(
     152             :       const double rest_mass_density) const override {
     153             :     return get(cold_eos_.specific_internal_energy_from_density(
     154             :                Scalar<double>{rest_mass_density})) +
     155             :            (min_temperature_) / (thermal_adiabatic_index_ - 1.0);
     156             :   }
     157             : 
     158             :   /// The upper bound of the specific internal energy that is valid for this EOS
     159             :   /// at the given rest mass density \f$\rho\f$
     160           1 :   double specific_internal_energy_upper_bound(
     161             :       const double /*rest_mass_density*/) const override {
     162             :     return std::numeric_limits<double>::max();
     163             :   }
     164             : 
     165             :   /// The lower bound of the specific enthalpy that is valid for this EOS
     166           1 :   double specific_enthalpy_lower_bound() const override {
     167             :     return cold_eos_.specific_enthalpy_lower_bound() +
     168             :            (thermal_adiabatic_index_ * min_temperature_) /
     169             :                (thermal_adiabatic_index_ - 1.0);
     170             :   }
     171             : 
     172             :   /// The lower bound of the temperature that is valid for this EOS.
     173             :   /// Non-zero lower bound could be set to impose floor on the specific
     174             :   /// internal energy.
     175           1 :   double temperature_lower_bound() const override { return min_temperature_; }
     176             : 
     177             :   /// The vacuum baryon mass for this EoS
     178           1 :   double baryon_mass() const override { return cold_eos_.baryon_mass(); }
     179             : 
     180             :  private:
     181             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(2)
     182             : 
     183           0 :   double specific_entropy_from_density_and_thermal_energy(
     184             :       double density, double energy) const;
     185             : 
     186           0 :   ColdEquationOfState cold_eos_;
     187           0 :   double thermal_adiabatic_index_ =
     188             :       std::numeric_limits<double>::signaling_NaN();
     189           0 :   double min_temperature_ = std::numeric_limits<double>::signaling_NaN();
     190             : };
     191             : 
     192             : /// \cond
     193             : template <typename ColdEquationOfState>
     194             : PUP::able::PUP_ID EquationsOfState::HybridEos<ColdEquationOfState>::my_PUP_ID =
     195             :     0;
     196             : /// \endcond
     197             : }  // namespace EquationsOfState