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             : #include <vector>
      16             : 
      17             : #include "DataStructures/Tensor/TypeAliases.hpp"
      18             : #include "Options/String.hpp"
      19             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"  // IWYU pragma: keep
      20             : #include "PointwiseFunctions/Hydro/Units.hpp"
      21             : #include "Utilities/Math.hpp"
      22             : #include "Utilities/Serialization/CharmPupable.hpp"
      23             : #include "Utilities/TMPL.hpp"
      24             : 
      25             : /// \cond
      26             : class DataVector;
      27             : 
      28             : /// \endcond
      29             : 
      30             : // IWYU pragma: no_forward_declare Tensor
      31             : 
      32             : namespace EquationsOfState {
      33             : 
      34             : /*!
      35             :  * \ingroup EquationsOfStateGroup
      36             :  * \brief An equation of state given by parametrized enthalpy
      37             :  *
      38             :  * This equation of state is determined as a function of \f$x =
      39             :  * \ln(\rho/\rho_0)\f$ where \f$\rho\f$ is the rest mass density and
      40             :  * \f$\rho_0\f$ is the provided reference density.
      41             :  * The pseudo-enthalpy \f$h \equiv (p + rho  + u)/rho\f$
      42             :  * is expanded as
      43             :  *
      44             :  * \f{equation}
      45             :  * h(x) = \sum_i a_i x^i + \sum_j b_j \sin(jkx) + c_j \cos(jkx)
      46             :  * \f}
      47             :  *
      48             :  * This form allows for convenient calculation of thermodynamic
      49             :  * quantities for a cold equation of state. For example
      50             :  *
      51             :  * \f{equation}
      52             :  * h(x) = \frac{d e} {d \rho} |_{x = \log(\rho/\rho_0)}
      53             :  * \f}
      54             :  *
      55             :  * where \f$e\f$ is the total energy density.  At the same time \f$ dx =
      56             :  * d\rho/\rho \f$ so \f$ \rho_0 e^x dx = d \rho \f$ Therefore,
      57             :  *
      58             :  * \f{equation}
      59             :  * e(x) - e(x_0) = \int_{x_0}^x h(x') e^{x'} dx '
      60             :  * \f}
      61             :  *
      62             :  * This can be computed analytically because
      63             :  *
      64             :  * \f{equation}
      65             :  *  \int a_i \frac{x^i}{i!} e^{x} dx = \sum_{j \leq i} a_i (-1)^{i-j}
      66             :  * \frac{(x)^{j}}{j!}
      67             :  * + C \f}
      68             :  *
      69             :  * and
      70             :  *
      71             :  * \f{equation}
      72             :  * \int b_j \sin(j k x) e^x dx = b_j e^x \frac{\sin(jkx) - j k \cos(jkx)}{j^2
      73             :  * k^2 + 1} \f}
      74             :  *
      75             :  * \f{equation}
      76             :  * \int c_j \cos(j k x) e^x dx = b_j e^x \frac{\cos(jkx) + j k \sin(jkx)}{j^2
      77             :  * k^2 + 1} \f}
      78             :  *
      79             :  * From this most other thermodynamic quantities can be computed
      80             :  * analytically
      81             :  *
      82             :  * The internal energy density
      83             :  * \f{equation}
      84             :  * \epsilon(x)\rho(x) = e(x)  - \rho(x)
      85             :  * \f}
      86             :  *
      87             :  * The pressure
      88             :  * \f{equation}
      89             :  * p(x) = \rho(x) h(x) - e(x)
      90             :  * \f}
      91             :  *
      92             :  * The derivative of the pressure with respect to the rest mass density
      93             :  * \f{equation}
      94             :  * \chi(x) = \frac{dp}{d\rho} |_{x = x(\rho)} = \frac{dh}{dx}
      95             :  * \f}
      96             :  *
      97             :  * Below the minimum density, a spectral parameterization
      98             :  * is used.
      99             :  *
     100             :  *
     101             :  *
     102             :  */
     103             : template <typename LowDensityEoS>
     104           1 : class Enthalpy : public EquationOfState<true, 1> {
     105             :  private:
     106           0 :   struct Coefficients {
     107           0 :     std::vector<double> polynomial_coefficients;
     108           0 :     std::vector<double> sin_coefficients;
     109           0 :     std::vector<double> cos_coefficients;
     110           0 :     double trig_scale;
     111           0 :     double reference_density;
     112           0 :     bool has_exponential_prefactor;
     113           0 :     double exponential_external_constant;
     114           0 :     Coefficients() = default;
     115           0 :     ~Coefficients() = default;
     116           0 :     Coefficients(const Coefficients& coefficients) = default;
     117           0 :     Coefficients(std::vector<double> in_polynomial_coefficients,
     118             :                  std::vector<double> in_sin_coefficients,
     119             :                  std::vector<double> in_cos_coefficients, double in_trig_scale,
     120             :                  double in_reference_density,
     121             :                  double in_exponential_constant =
     122             :                      std::numeric_limits<double>::quiet_NaN());
     123           0 :     bool operator==(const Coefficients& rhs) const;
     124             : 
     125           0 :     Enthalpy<LowDensityEoS>::Coefficients compute_exponential_integral(
     126             :         const std::pair<double, double>& initial_condition,
     127             :         const double minimum_density);
     128           0 :     Enthalpy<LowDensityEoS>::Coefficients compute_derivative();
     129           0 :     void pup(PUP::er& p);
     130             :   };
     131             : 
     132             :  public:
     133           0 :   static constexpr size_t thermodynamic_dim = 1;
     134           0 :   static constexpr bool is_relativistic = true;
     135             : 
     136           0 :   static std::string name() {
     137             :     return "Enthalpy(" + pretty_type::name<LowDensityEoS>() + ")";
     138             :   }
     139             : 
     140           0 :   struct ReferenceDensity {
     141           0 :     using type = double;
     142           0 :     static constexpr Options::String help = {"Reference density rho_0"};
     143           0 :     static double lower_bound() { return 0.0; }
     144             :   };
     145             : 
     146           0 :   struct MinimumDensity {
     147           0 :     using type = double;
     148           0 :     static constexpr Options::String help = {
     149             :         "Minimum valid density rho_min,"
     150             :         " for this parametrization"};
     151           0 :     static double lower_bound() { return 0.0; }
     152             :   };
     153           0 :   struct MaximumDensity {
     154           0 :     using type = double;
     155           0 :     static constexpr Options::String help = {"Maximum density for this EoS"};
     156           0 :     static double lower_bound() { return 0.0; }
     157             :   };
     158             : 
     159           0 :   struct PolynomialCoefficients {
     160           0 :     using type = std::vector<double>;
     161           0 :     static constexpr Options::String help = {"Polynomial coefficients a_i"};
     162             :   };
     163             : 
     164           0 :   struct TrigScaling {
     165           0 :     using type = double;
     166           0 :     static constexpr Options::String help = {
     167             :         "Fundamental wavenumber of trig "
     168             :         "functions, k"};
     169           0 :     static double lower_bound() { return 0.0; }
     170             :   };
     171             : 
     172           0 :   struct SinCoefficients {
     173           0 :     using type = std::vector<double>;
     174           0 :     static constexpr Options::String help = {"Sine coefficients b_j"};
     175             :   };
     176           0 :   struct CosCoefficients {
     177           0 :     using type = std::vector<double>;
     178           0 :     static constexpr Options::String help = {"Cosine coefficients c_j"};
     179             :   };
     180           0 :   struct StitchedLowDensityEoS {
     181           0 :     using type = LowDensityEoS;
     182           0 :     static std::string name() {
     183             :       return pretty_type::short_name<LowDensityEoS>();
     184             :     }
     185           0 :     static constexpr Options::String help = {
     186             :         "Low density EoS stitched at the MinimumDensity"};
     187             :   };
     188             : 
     189           0 :   struct TransitionDeltaEpsilon {
     190           0 :     using type = double;
     191           0 :     static constexpr Options::String help = {
     192             :         "the change in internal energy across the low-"
     193             :         "to-high-density transition, generically 0.0"};
     194           0 :     static double lower_bound() { return 0.0; }
     195             :   };
     196             : 
     197           0 :   static constexpr Options::String help = {
     198             :       "An EoS with a parametrized value h(log(rho/rho_0)) with h the specific "
     199             :       "enthalpy and rho the baryon rest mass density.  The enthalpy is "
     200             :       "expanded as a sum of polynomial terms and trigonometric corrections. "
     201             :       "let x = log(rho/rho_0) in"
     202             :       "h(x) = \\sum_i a_ix^i + \\sum_j b_jsin(k * j * x) + c_jcos(k * j * x) "
     203             :       "Note that rho(x)(1+epsilon(x)) = int_0^x e^x' h((x') dx' can be "
     204             :       "computed "
     205             :       "analytically, and therefore so can "
     206             :       "P(x) = rho(x) * (h(x) - (1 + epsilon(x))) "};
     207             : 
     208           0 :   using options =
     209             :       tmpl::list<ReferenceDensity, MaximumDensity, MinimumDensity, TrigScaling,
     210             :                  PolynomialCoefficients, SinCoefficients, CosCoefficients,
     211             :                  StitchedLowDensityEoS, TransitionDeltaEpsilon>;
     212             : 
     213           0 :   Enthalpy() = default;
     214           0 :   Enthalpy(const Enthalpy&) = default;
     215           0 :   Enthalpy& operator=(const Enthalpy&) = default;
     216           0 :   Enthalpy(Enthalpy&&) = default;
     217           0 :   Enthalpy& operator=(Enthalpy&&) = default;
     218           0 :   ~Enthalpy() override = default;
     219             : 
     220           0 :   Enthalpy(double reference_density, double max_density, double min_density,
     221             :            double trig_scale,
     222             :            const std::vector<double>& polynomial_coefficients,
     223             :            const std::vector<double>& sin_coefficients,
     224             :            const std::vector<double>& cos_coefficients,
     225             :            const LowDensityEoS& low_density_eos,
     226             :            const double transition_delta_epsilon);
     227             : 
     228           0 :   std::unique_ptr<EquationOfState<true, 1>> get_clone() const override;
     229             : 
     230           0 :   std::unique_ptr<EquationOfState<true, 3>> promote_to_3d_eos() const override;
     231             : 
     232           0 :   std::unique_ptr<EquationOfState<true, 2>> promote_to_2d_eos() const override;
     233             : 
     234           0 :   bool is_equal(const EquationOfState<true, 1>& rhs) const override;
     235             : 
     236           0 :   bool operator==(const Enthalpy<LowDensityEoS>& rhs) const;
     237             : 
     238           0 :   bool operator!=(const Enthalpy<LowDensityEoS>& rhs) const;
     239             : 
     240             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(Enthalpy, 1)
     241             : 
     242           0 :   WRAPPED_PUPable_decl_base_template(  // NOLINT
     243             :       SINGLE_ARG(EquationOfState<true, 1>), Enthalpy);
     244             : 
     245             :   /// The lower bound of the rest mass density that is valid for this EOS
     246           1 :   double rest_mass_density_lower_bound() const override { return 0.0; }
     247             : 
     248             :   /// The upper bound of the rest mass density that is valid for this EOS
     249           1 :   double rest_mass_density_upper_bound() const override {
     250             :     return std::numeric_limits<double>::max();
     251             :   }
     252             : 
     253             :   /// The lower bound of the specific internal energy that is valid for this EOS
     254             :   /// at the given rest mass density \f$\rho\f$
     255           1 :   double specific_internal_energy_lower_bound(
     256             :       const double /* rest_mass_density */) const override {
     257             :     return 0.0;
     258             :   }
     259             : 
     260             :   /// The upper bound of the specific internal energy that is valid for this EOS
     261             :   /// at the given rest mass density \f$\rho\f$
     262           1 :   double specific_internal_energy_upper_bound(
     263             :       const double /* rest_mass_density */) const override {
     264             :     return std::numeric_limits<double>::max();
     265             :   }
     266             : 
     267             :   /// The lower bound of the specific enthalpy that is valid for this EOS
     268           1 :   double specific_enthalpy_lower_bound() const override { return 1.0; }
     269             : 
     270             :   /// The vacuum baryon mass for this EoS
     271           1 :   double baryon_mass() const override { return low_density_eos_.baryon_mass(); }
     272             : 
     273             :  private:
     274             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(1)
     275             : 
     276             :   SPECTRE_ALWAYS_INLINE
     277           0 :   bool in_low_density_domain(const double rest_mass_density) const {
     278             :     return rest_mass_density < minimum_density_;
     279             :   }
     280             : 
     281           0 :   double x_from_density(const double rest_mass_density) const;
     282           0 :   double density_from_x(const double x) const;
     283           0 :   double energy_density_from_log_density(const double x,
     284             :                                          const double rest_mass_density) const;
     285           0 :   static double evaluate_coefficients(
     286             :       const Enthalpy::Coefficients& coefficients, const double x,
     287             :       const double exponential_prefactor =
     288             :           std::numeric_limits<double>::signaling_NaN());
     289           0 :   static Enthalpy::Coefficients compute_pressure_coefficients(
     290             :       const typename Enthalpy::Coefficients& enthalpy,
     291             :       const typename Enthalpy::Coefficients& energy_density);
     292             : 
     293           0 :   double chi_from_density(const double rest_mass_density) const;
     294           0 :   double specific_internal_energy_from_density(
     295             :       const double rest_mass_density) const;
     296           0 :   double specific_enthalpy_from_density(const double rest_mass_density) const;
     297           0 :   double pressure_from_density(const double rest_mass_density) const;
     298           0 :   double pressure_from_log_density(const double x,
     299             :                                    const double rest_mass_density) const;
     300           0 :   double rest_mass_density_from_enthalpy(const double specific_enthalpy) const;
     301             : 
     302           0 :   double reference_density_ = std::numeric_limits<double>::signaling_NaN();
     303           0 :   double minimum_density_ = std::numeric_limits<double>::signaling_NaN();
     304           0 :   double maximum_density_ = std::numeric_limits<double>::signaling_NaN();
     305           0 :   double minimum_enthalpy_ = std::numeric_limits<double>::signaling_NaN();
     306             : 
     307           0 :   LowDensityEoS low_density_eos_;
     308           0 :   Coefficients coefficients_;
     309           0 :   Coefficients exponential_integral_coefficients_;
     310           0 :   Coefficients derivative_coefficients_;
     311           0 :   Coefficients pressure_coefficients_;
     312             : };
     313             : 
     314             : }  // namespace EquationsOfState