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"  // IWYU pragma: keep
      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             : // IWYU pragma: no_forward_declare Tensor
      28             : 
      29             : namespace EquationsOfState {
      30             : /*!
      31             :  * \ingroup EquationsOfStateGroup
      32             :  * \brief Equation of state for a piecewise polytropic fluid
      33             :  *
      34             :  * A piecewise polytropic equation of state \f$p=K_i\rho^{\Gamma_i}\f$ where
      35             :  *  \f$K_i\f$ is the polytropic constant and \f$\Gamma_i\f$ is the polytropic
      36             :  * exponent. Here the subscript \f$i\f$ indicates two pairs of constants and
      37             :  *  exponents which characterize `the stiffness' of the matter at low and high
      38             :  *  densities.  For a given density, the polytropic exponent is related to the
      39             :  *  polytropic index \f$N_p\f$ by \f$N_p=1/(\Gamma-1)\f$.  For posterity,
      40             :  *  this two piece polytrope has been used in toy models of CCSNe (e.g.,
      41             :  *  \cite Dimmelmeier2001 ) and could be extended to a general "M" number of
      42             :  * parts for simplified equations of state for neutron stars (e.g.,
      43             :  *  \cite OBoyle2020 ). For a reference to a general piecewise polytrope, see
      44             :  * Section 2.4.7 of \cite RezzollaBook.
      45             :  */
      46             : template <bool IsRelativistic>
      47           1 : class PiecewisePolytropicFluid : public EquationOfState<IsRelativistic, 1> {
      48             :  public:
      49           0 :   static constexpr size_t thermodynamic_dim = 1;
      50           0 :   static constexpr bool is_relativistic = IsRelativistic;
      51             : 
      52             :   /// The density demarcating the high and low density descriptions of the
      53             :   /// fluid.
      54           1 :   struct PiecewisePolytropicTransitionDensity {
      55           0 :     using type = double;
      56           0 :     static constexpr Options::String help = {
      57             :         "Density below (above) which, the matter is described by a low (high) "
      58             :         "density polytropic fluid."};
      59           0 :     static double lower_bound() { return 0.0; }
      60             :   };
      61             : 
      62             :   /// The constant \f$K\f$ scaling the low density material
      63             :   /// \f$p=K\rho^{\Gamma}\f$.
      64             :   ///
      65             :   /// Note, by enforcing pressure continuity at the transition density
      66             :   /// \f$\bar{\rho}\f$, the high density constant \f$K_{high}\f$ is given
      67             :   /// as \f$K_{high} = K_{low} (\bar{\rho})^{\Gamma_{low} - \Gamma_{high}}\f$.
      68           1 :   struct PolytropicConstantLow {
      69           0 :     using type = double;
      70           0 :     static constexpr Options::String help = {
      71             :         "Polytropic constant K for lower"
      72             :         " density material"};
      73           0 :     static double lower_bound() { return 0.0; }
      74             :   };
      75             : 
      76             :   /// The exponent \f$\Gamma\f$, scaling the low density material
      77             :   /// \f$p=K\rho^{\Gamma}\f$.
      78           1 :   struct PolytropicExponentLow {
      79           0 :     using type = double;
      80           0 :     static constexpr Options::String help = {
      81             :         "Polytropic exponent for lower"
      82             :         " density material."};
      83           0 :     static double lower_bound() { return 1.0; }
      84             :   };
      85             : 
      86             :   /// The exponent \f$\Gamma\f$, scaling the high density material
      87             :   /// \f$p=K\rho^{\Gamma}\f$.
      88           1 :   struct PolytropicExponentHigh {
      89           0 :     using type = double;
      90           0 :     static constexpr Options::String help = {
      91             :         "Polytropic exponent for higher"
      92             :         " density material."};
      93           0 :     static double lower_bound() { return 1.0; }
      94             :   };
      95             : 
      96           0 :   static constexpr Options::String help = {
      97             :       "A piecewise polytropic fluid equation of state.\n"
      98             :       "The pressure is related to the rest mass density by p = K_i rho ^ "
      99             :       "Gamma_i, "
     100             :       "where p is the pressure, rho is the rest mass density, K_i is the "
     101             :       "polytropic constant either describing the low or high density material, "
     102             :       "and Gamma_i is the polytropic exponent for the low or high density "
     103             :       "material. The polytropic index N_i is defined as Gamma_i = 1 + 1 / N_i."
     104             :       "  The subscript `i' refers to different pairs of Gamma and K that can"
     105             :       " describe either low or high density material."};
     106             : 
     107           0 :   using options =
     108             :       tmpl::list<PiecewisePolytropicTransitionDensity, PolytropicConstantLow,
     109             :                  PolytropicExponentLow, PolytropicExponentHigh>;
     110             : 
     111           0 :   PiecewisePolytropicFluid() = default;
     112           0 :   PiecewisePolytropicFluid(const PiecewisePolytropicFluid&) = default;
     113           0 :   PiecewisePolytropicFluid& operator=(const PiecewisePolytropicFluid&) =
     114             :       default;
     115           0 :   PiecewisePolytropicFluid(PiecewisePolytropicFluid&&) = default;
     116           0 :   PiecewisePolytropicFluid& operator=(PiecewisePolytropicFluid&&) = default;
     117           0 :   ~PiecewisePolytropicFluid() override = default;
     118             : 
     119           0 :   PiecewisePolytropicFluid(double transition_density,
     120             :                            double polytropic_constant_lo,
     121             :                            double polytropic_exponent_lo,
     122             :                            double polytropic_exponent_hi);
     123             : 
     124           0 :   std::unique_ptr<EquationOfState<IsRelativistic, 1>> get_clone()
     125             :       const override;
     126             : 
     127           1 :   std::unique_ptr<EquationOfState<IsRelativistic, 3>> promote_to_3d_eos()
     128             :       const override;
     129             : 
     130           1 :   std::unique_ptr<EquationOfState<IsRelativistic, 2>> promote_to_2d_eos()
     131             :       const override;
     132             : 
     133           0 :   bool is_equal(const EquationOfState<IsRelativistic, 1>& rhs) const override;
     134             : 
     135           0 :   bool operator==(const PiecewisePolytropicFluid<IsRelativistic>& rhs) const;
     136             : 
     137           0 :   bool operator!=(const PiecewisePolytropicFluid<IsRelativistic>& rhs) const;
     138             : 
     139             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(PiecewisePolytropicFluid, 1)
     140             : 
     141           0 :   WRAPPED_PUPable_decl_base_template(  // NOLINT
     142             :       SINGLE_ARG(EquationOfState<IsRelativistic, 1>), PiecewisePolytropicFluid);
     143             : 
     144             :   /// The lower bound of the rest mass density that is valid for this EOS
     145           1 :   double rest_mass_density_lower_bound() const override { return 0.0; }
     146             : 
     147             :   /// The upper bound of the rest mass density that is valid for this EOS
     148           1 :   double rest_mass_density_upper_bound() const override;
     149             : 
     150             :   /// The lower bound of the specific internal energy that is valid for this EOS
     151             :   /// at the given rest mass density \f$\rho\f$
     152           1 :   double specific_internal_energy_lower_bound(
     153             :       const double /* rest_mass_density */) const override {
     154             :     return 0.0;
     155             :   }
     156             : 
     157             :   /// The upper bound of the specific internal energy that is valid for this EOS
     158             :   /// at the given rest mass density \f$\rho\f$
     159           1 :   double specific_internal_energy_upper_bound(
     160             :       const double /* rest_mass_density */) const override {
     161             :     return std::numeric_limits<double>::max();
     162             :   }
     163             : 
     164             :   /// The lower bound of the specific enthalpy that is valid for this EOS
     165           1 :   double specific_enthalpy_lower_bound() const override {
     166             :     return IsRelativistic ? 1.0 : 0.0;
     167             :   }
     168             : 
     169             :   /// The vacuum baryon mass for this EoS
     170           1 :   double baryon_mass() const override {
     171             :     return hydro::units::geometric::default_baryon_mass;
     172             :   }
     173             : 
     174             :  private:
     175             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(1)
     176             : 
     177           0 :   double transition_density_ = std::numeric_limits<double>::signaling_NaN();
     178           0 :   double transition_pressure_ = std::numeric_limits<double>::signaling_NaN();
     179           0 :   double transition_spec_eint_ = std::numeric_limits<double>::signaling_NaN();
     180           0 :   double polytropic_constant_lo_ = std::numeric_limits<double>::signaling_NaN();
     181           0 :   double polytropic_exponent_lo_ = std::numeric_limits<double>::signaling_NaN();
     182           0 :   double polytropic_constant_hi_ = std::numeric_limits<double>::signaling_NaN();
     183           0 :   double polytropic_exponent_hi_ = std::numeric_limits<double>::signaling_NaN();
     184             : };
     185             : 
     186             : /// \cond
     187             : template <bool IsRelativistic>
     188             : PUP::able::PUP_ID
     189             :     EquationsOfState::PiecewisePolytropicFluid<IsRelativistic>::my_PUP_ID = 0;
     190             : /// \endcond
     191             : }  // namespace EquationsOfState