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