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