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