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