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 <limits>
      14             : #include <pup.h>
      15             : 
      16             : #include "DataStructures/Tensor/TypeAliases.hpp"
      17             : #include "Options/String.hpp"
      18             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      19             : #include "PointwiseFunctions/Hydro/Units.hpp"
      20             : #include "Utilities/Serialization/CharmPupable.hpp"
      21             : #include "Utilities/TMPL.hpp"
      22             : 
      23             : /// \cond
      24             : class DataVector;
      25             : /// \endcond
      26             : 
      27             : namespace EquationsOfState {
      28             : /*!
      29             :  * \ingroup EquationsOfStateGroup
      30             :  * \brief Equation of state for a polytropic fluid
      31             :  *
      32             :  * A polytropic equation of state \f$p=K\rho^{\Gamma}\f$ where \f$K\f$ is the
      33             :  * polytropic constant and \f$\Gamma\f$ is the polytropic exponent. The
      34             :  * polytropic exponent is related to the polytropic index \f$N_p\f$ by
      35             :  * \f$N_p=1/(\Gamma-1)\f$.
      36             :  *
      37             :  * We also have
      38             :  *
      39             :  * \f{align}{
      40             :  * \epsilon&=\frac{K\rho^{\Gamma-1}}{\Gamma-1}\\
      41             :  * h&=1+\epsilon+\frac{p}{\rho}=1+\frac{K\Gamma}{\Gamma-1}\rho^{\Gamma-1} \\
      42             :  * T&=0 \\
      43             :  * c_s^2&=\frac{\Gamma p}{\rho h}
      44             :  * =\frac{\Gamma(\Gamma-1)p}{\rho(\Gamma-1)+\Gamma p}
      45             :  * =\left(\frac{1}{\Gamma K\rho^{\Gamma-1}}+\frac{1}{\Gamma-1}\right)^{-1}
      46             :  * \f}
      47             :  */
      48             : template <bool IsRelativistic>
      49           1 : class PolytropicFluid : public EquationOfState<IsRelativistic, 1> {
      50             :  public:
      51           0 :   static constexpr size_t thermodynamic_dim = 1;
      52           0 :   static constexpr bool is_relativistic = IsRelativistic;
      53             : 
      54           0 :   struct PolytropicConstant {
      55           0 :     using type = double;
      56           0 :     static constexpr Options::String help = {"Polytropic constant K"};
      57           0 :     static double lower_bound() { return 0.0; }
      58             :   };
      59             : 
      60           0 :   struct PolytropicExponent {
      61           0 :     using type = double;
      62           0 :     static constexpr Options::String help = {"Polytropic exponent Gamma"};
      63           0 :     static double lower_bound() { return 1.0; }
      64             :   };
      65             : 
      66           0 :   static constexpr Options::String help = {
      67             :       "A polytropic fluid equation of state.\n"
      68             :       "The pressure is related to the rest mass density by p = K rho ^ Gamma, "
      69             :       "where p is the pressure, rho is the rest mass density, K is the "
      70             :       "polytropic constant, and Gamma is the polytropic exponent. The "
      71             :       "polytropic index N is defined as Gamma = 1 + 1 / N."};
      72             : 
      73           0 :   using options = tmpl::list<PolytropicConstant, PolytropicExponent>;
      74             : 
      75           0 :   PolytropicFluid() = default;
      76           0 :   PolytropicFluid(const PolytropicFluid&) = default;
      77           0 :   PolytropicFluid& operator=(const PolytropicFluid&) = default;
      78           0 :   PolytropicFluid(PolytropicFluid&&) = default;
      79           0 :   PolytropicFluid& operator=(PolytropicFluid&&) = default;
      80           0 :   ~PolytropicFluid() override = default;
      81             : 
      82           0 :   PolytropicFluid(double polytropic_constant, double polytropic_exponent);
      83             : 
      84           0 :   std::unique_ptr<EquationOfState<IsRelativistic, 1>> get_clone()
      85             :       const override;
      86             : 
      87           1 :   std::unique_ptr<EquationOfState<IsRelativistic, 3>> promote_to_3d_eos()
      88             :       const override;
      89             : 
      90           1 :   std::unique_ptr<EquationOfState<IsRelativistic, 2>> promote_to_2d_eos()
      91             :       const override;
      92             : 
      93           0 :   bool is_equal(const EquationOfState<IsRelativistic, 1>& rhs) const override;
      94             : 
      95           0 :   bool operator==(const PolytropicFluid<IsRelativistic>& rhs) const;
      96             : 
      97           0 :   bool operator!=(const PolytropicFluid<IsRelativistic>& rhs) const;
      98             : 
      99             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(PolytropicFluid, 1)
     100             : 
     101           0 :   WRAPPED_PUPable_decl_base_template(  // NOLINT
     102             :       SINGLE_ARG(EquationOfState<IsRelativistic, 1>), PolytropicFluid);
     103             : 
     104             :   /// The lower bound of the rest mass density that is valid for this EOS
     105           1 :   double rest_mass_density_lower_bound() const override { return 0.0; }
     106             : 
     107             :   /// The upper bound of the rest mass density that is valid for this EOS
     108           1 :   double rest_mass_density_upper_bound() const override;
     109             : 
     110             :   /// The lower bound of the specific enthalpy that is valid for this EOS
     111           1 :   double specific_enthalpy_lower_bound() const override {
     112             :     return IsRelativistic ? 1.0 : 0.0;
     113             :   }
     114             : 
     115             :   /// The lower bound of the specific internal energy that is valid for this EOS
     116           1 :   double specific_internal_energy_lower_bound() const override { return 0.0; }
     117             : 
     118             :   /// The upper bound of the specific internal energy that is valid for this EOS
     119           1 :   double specific_internal_energy_upper_bound() const override;
     120             : 
     121             :   /// The vacuum baryon mass for this EoS
     122           1 :   double baryon_mass() const override {
     123             :     return hydro::units::geometric::default_baryon_mass;
     124             :   }
     125             : 
     126             :  private:
     127             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(1)
     128             : 
     129           0 :   double polytropic_constant_ = std::numeric_limits<double>::signaling_NaN();
     130           0 :   double polytropic_exponent_ = std::numeric_limits<double>::signaling_NaN();
     131             : };
     132             : 
     133             : /// \cond
     134             : template <bool IsRelativistic>
     135             : PUP::able::PUP_ID EquationsOfState::PolytropicFluid<IsRelativistic>::my_PUP_ID =
     136             :     0;
     137             : /// \endcond
     138             : }  // namespace EquationsOfState