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             : #include <vector>
      17             : 
      18             : #include "DataStructures/Tensor/TypeAliases.hpp"
      19             : #include "Options/String.hpp"
      20             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"  // IWYU pragma: keep
      21             : #include "PointwiseFunctions/Hydro/Units.hpp"
      22             : #include "Utilities/Math.hpp"
      23             : #include "Utilities/Serialization/CharmPupable.hpp"
      24             : #include "Utilities/TMPL.hpp"
      25             : 
      26             : /// \cond
      27             : class DataVector;
      28             : /// \endcond
      29             : 
      30             : // IWYU pragma: no_forward_declare Tensor
      31             : 
      32             : namespace EquationsOfState {
      33             : /*!
      34             :  * \ingroup EquationsOfStateGroup
      35             :  * \brief A spectral equation of state
      36             :  *
      37             :  * This equation of state is determined as a function of \f$x =
      38             :  * \ln(\rho/\rho_0)\f$ where \f$\rho\f$ is the rest mass density and
      39             :  * \f$\rho_0\f$ is the provided reference density.  The adiabatic
      40             :  * index \f$\Gamma(x)\f$ is defined such that
      41             :  * \f{equation}{
      42             :  * \frac{d \ln p}{dx} = \Gamma(x) = \sum_{n=0}^N
      43             :  * \gamma_n x^n
      44             :  * \f}
      45             :  *
      46             :  * for the set of spectral coefficinets \f$\gamma_n\f$ when
      47             :  * \f$0 < x < x_u = \ln(\rho_u/\rho_0)\f$, where \f$\rho_u\f$ is the provided
      48             :  * upper density.
      49             :  *
      50             :  * For \f$ x < 0 \f$, \f$ \Gamma(x) = \gamma_0 \f$.
      51             :  *
      52             :  * For \f$ x > x_u \f$, \f$ \Gamma(x) = \Gamma(x_u) \f$
      53             :  *
      54             :  *
      55             :  */
      56           1 : class Spectral : public EquationOfState<true, 1> {
      57             :  public:
      58           0 :   static constexpr size_t thermodynamic_dim = 1;
      59           0 :   static constexpr bool is_relativistic = true;
      60             : 
      61           0 :   struct ReferenceDensity {
      62           0 :     using type = double;
      63           0 :     static constexpr Options::String help = {"Reference density rho_0"};
      64           0 :     static double lower_bound() { return 0.0; }
      65             :   };
      66             : 
      67           0 :   struct ReferencePressure {
      68           0 :     using type = double;
      69           0 :     static constexpr Options::String help = {"Reference pressure p_0"};
      70           0 :     static double lower_bound() { return 0.0; }
      71             :   };
      72             : 
      73           0 :   struct Coefficients {
      74           0 :     using type = std::vector<double>;
      75           0 :     static constexpr Options::String help = {"Spectral coefficients gamma_i"};
      76             :   };
      77             : 
      78           0 :   struct UpperDensity {
      79           0 :     using type = double;
      80           0 :     static constexpr Options::String help = {"Upper density rho_u"};
      81           0 :     static double lower_bound() { return 0.0; }
      82             :   };
      83             : 
      84           0 :   static constexpr Options::String help = {
      85             :       "A spectral equation of state.  Defining x = log(rho/rho_0), Gamma(x) = "
      86             :       "Sum_i gamma_i x^i, then the pressure is determined from d(log P)/dx = "
      87             :       "Gamma(x) for x > 0.  For x < 0 the EOS is a polytrope with "
      88             :       "Gamma(x)=Gamma(0).  For x > x_u = log(rho_u/rho_0), Gamma(x) = "
      89             :       "Gamma(x_u).\n"
      90             :       "To get smooth equations of state, it is recommended that the second "
      91             :       "and third supplied coefficient should be 0. It is up to the user to "
      92             :       "choose coefficients that are physically reasonable, e.g. that "
      93             :       "satisfy causality."};
      94             : 
      95           0 :   using options = tmpl::list<ReferenceDensity, ReferencePressure, Coefficients,
      96             :                              UpperDensity>;
      97             : 
      98           0 :   Spectral() = default;
      99           0 :   Spectral(const Spectral&) = default;
     100           0 :   Spectral& operator=(const Spectral&) = default;
     101           0 :   Spectral(Spectral&&) = default;
     102           0 :   Spectral& operator=(Spectral&&) = default;
     103           0 :   ~Spectral() override = default;
     104             : 
     105           0 :   Spectral(double reference_density, double reference_pressure,
     106             :            std::vector<double> coefficients, double upper_density);
     107             : 
     108             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(Spectral, 1)
     109             : 
     110           0 :   std::unique_ptr<EquationOfState<true, 1>> get_clone() const override;
     111             : 
     112           0 :   std::unique_ptr<EquationOfState<true, 3>> promote_to_3d_eos() const override;
     113             : 
     114           0 :   std::unique_ptr<EquationOfState<true, 2>> promote_to_2d_eos() const override;
     115             : 
     116           0 :   bool operator==(const Spectral& rhs) const;
     117             : 
     118           0 :   bool operator!=(const Spectral& rhs) const;
     119             : 
     120           0 :   bool is_equal(const EquationOfState<true, 1>& rhs) const override;
     121             : 
     122           0 :   WRAPPED_PUPable_decl_base_template(  // NOLINT
     123             :       SINGLE_ARG(EquationOfState<true, 1>), Spectral);
     124             : 
     125             :   /// The lower bound of the rest mass density that is valid for this EOS
     126           1 :   double rest_mass_density_lower_bound() const override { return 0.0; }
     127             : 
     128             :   /// The upper bound of the rest mass density that is valid for this EOS
     129           1 :   double rest_mass_density_upper_bound() const override {
     130             :     return std::numeric_limits<double>::max();
     131             :   }
     132             : 
     133             :   /// The lower bound of the specific internal energy that is valid for this EOS
     134             :   /// at the given rest mass density \f$\rho\f$
     135           1 :   double specific_internal_energy_lower_bound(
     136             :       const double /* rest_mass_density */) const override {
     137             :     return 0.0;
     138             :   }
     139             : 
     140             :   /// The upper bound of the specific internal energy that is valid for this EOS
     141             :   /// at the given rest mass density \f$\rho\f$
     142           1 :   double specific_internal_energy_upper_bound(
     143             :       const double /* rest_mass_density */) const override {
     144             :     return std::numeric_limits<double>::max();
     145             :   }
     146             : 
     147             :   /// The lower bound of the specific enthalpy that is valid for this EOS
     148           1 :   double specific_enthalpy_lower_bound() const override { return 1.0; }
     149             : 
     150             :   /// The vacuum baryon mass for this EoS
     151           1 :   double baryon_mass() const override {
     152             :     return hydro::units::geometric::default_baryon_mass;
     153             :   }
     154             : 
     155             :  private:
     156             :   EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(1)
     157             : 
     158           0 :   double gamma(const double x) const;
     159           0 :   double integral_of_gamma(const double x) const;
     160           0 :   double chi_from_density(const double density) const;
     161           0 :   double specific_internal_energy_from_density(const double density) const;
     162           0 :   double specific_enthalpy_from_density(const double density) const;
     163           0 :   double pressure_from_density(const double density) const;
     164           0 :   double pressure_from_log_density(const double x) const;
     165           0 :   double rest_mass_density_from_enthalpy(const double specific_enthalpy) const;
     166             : 
     167           0 :   double reference_density_ = std::numeric_limits<double>::signaling_NaN();
     168           0 :   double reference_pressure_ = std::numeric_limits<double>::signaling_NaN();
     169           0 :   std::vector<double> integral_coefficients_{};
     170           0 :   std::vector<double> gamma_coefficients_{};
     171           0 :   double x_max_ = std::numeric_limits<double>::signaling_NaN();
     172           0 :   double gamma_of_x_max_ = std::numeric_limits<double>::signaling_NaN();
     173           0 :   double integral_of_gamma_of_x_max_ =
     174             :       std::numeric_limits<double>::signaling_NaN();
     175           0 :   std::vector<double> table_of_specific_energies_{};
     176             :   // Information for Gaussian quadrature
     177           0 :   size_t number_of_quadrature_coefs_ =
     178             :       std::numeric_limits<size_t>::signaling_NaN();
     179           0 :   std::vector<double> quadrature_weights_{};
     180           0 :   std::vector<double> quadrature_points_{};
     181             : };
     182             : 
     183             : }  // namespace EquationsOfState