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
|
