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 <limits>
14 : #include <pup.h>
15 : #include <vector>
16 :
17 : #include "DataStructures/Tensor/TypeAliases.hpp"
18 : #include "Options/String.hpp"
19 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
20 : #include "PointwiseFunctions/Hydro/Units.hpp"
21 : #include "Utilities/Math.hpp"
22 : #include "Utilities/Serialization/CharmPupable.hpp"
23 : #include "Utilities/TMPL.hpp"
24 :
25 : /// \cond
26 : class DataVector;
27 :
28 : /// \endcond
29 :
30 : namespace EquationsOfState {
31 :
32 : /*!
33 : * \ingroup EquationsOfStateGroup
34 : * \brief An equation of state given by parametrized enthalpy
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.
39 : * The pseudo-enthalpy \f$h \equiv (p + rho + u)/rho\f$
40 : * is expanded as
41 : *
42 : * \f{equation}
43 : * h(x) = \sum_i a_i x^i + \sum_j b_j \sin(jkx) + c_j \cos(jkx)
44 : * \f}
45 : *
46 : * This form allows for convenient calculation of thermodynamic
47 : * quantities for a cold equation of state. For example
48 : *
49 : * \f{equation}
50 : * h(x) = \frac{d e} {d \rho} |_{x = \log(\rho/\rho_0)}
51 : * \f}
52 : *
53 : * where \f$e\f$ is the total energy density. At the same time \f$ dx =
54 : * d\rho/\rho \f$ so \f$ \rho_0 e^x dx = d \rho \f$ Therefore,
55 : *
56 : * \f{equation}
57 : * e(x) - e(x_0) = \int_{x_0}^x h(x') e^{x'} dx '
58 : * \f}
59 : *
60 : * This can be computed analytically because
61 : *
62 : * \f{equation}
63 : * \int a_i \frac{x^i}{i!} e^{x} dx = \sum_{j \leq i} a_i (-1)^{i-j}
64 : * \frac{(x)^{j}}{j!}
65 : * + C \f}
66 : *
67 : * and
68 : *
69 : * \f{equation}
70 : * \int b_j \sin(j k x) e^x dx = b_j e^x \frac{\sin(jkx) - j k \cos(jkx)}{j^2
71 : * k^2 + 1} \f}
72 : *
73 : * \f{equation}
74 : * \int c_j \cos(j k x) e^x dx = b_j e^x \frac{\cos(jkx) + j k \sin(jkx)}{j^2
75 : * k^2 + 1} \f}
76 : *
77 : * From this most other thermodynamic quantities can be computed
78 : * analytically
79 : *
80 : * The internal energy density
81 : * \f{equation}
82 : * \epsilon(x)\rho(x) = e(x) - \rho(x)
83 : * \f}
84 : *
85 : * The pressure
86 : * \f{equation}
87 : * p(x) = \rho(x) h(x) - e(x)
88 : * \f}
89 : *
90 : * The derivative of the pressure with respect to the rest mass density
91 : * \f{equation}
92 : * \chi(x) = \frac{dp}{d\rho} |_{x = x(\rho)} = \frac{dh}{dx}
93 : * \f}
94 : *
95 : * Below the minimum density, a spectral parameterization
96 : * is used.
97 : *
98 : *
99 : *
100 : */
101 : template <typename LowDensityEoS>
102 1 : class Enthalpy : public EquationOfState<true, 1> {
103 : private:
104 0 : struct Coefficients {
105 0 : std::vector<double> polynomial_coefficients;
106 0 : std::vector<double> sin_coefficients;
107 0 : std::vector<double> cos_coefficients;
108 0 : double trig_scale;
109 0 : double reference_density;
110 0 : bool has_exponential_prefactor;
111 0 : double exponential_external_constant;
112 0 : Coefficients() = default;
113 0 : ~Coefficients() = default;
114 0 : Coefficients(const Coefficients& coefficients) = default;
115 0 : Coefficients(std::vector<double> in_polynomial_coefficients,
116 : std::vector<double> in_sin_coefficients,
117 : std::vector<double> in_cos_coefficients, double in_trig_scale,
118 : double in_reference_density,
119 : double in_exponential_constant =
120 : std::numeric_limits<double>::quiet_NaN());
121 0 : bool operator==(const Coefficients& rhs) const;
122 :
123 0 : Enthalpy<LowDensityEoS>::Coefficients compute_exponential_integral(
124 : const std::pair<double, double>& initial_condition,
125 : const double minimum_density);
126 0 : Enthalpy<LowDensityEoS>::Coefficients compute_derivative();
127 0 : void pup(PUP::er& p);
128 : };
129 :
130 : public:
131 0 : static constexpr size_t thermodynamic_dim = 1;
132 0 : static constexpr bool is_relativistic = true;
133 :
134 0 : static std::string name() {
135 : return "Enthalpy(" + pretty_type::name<LowDensityEoS>() + ")";
136 : }
137 :
138 0 : struct ReferenceDensity {
139 0 : using type = double;
140 0 : static constexpr Options::String help = {"Reference density rho_0"};
141 0 : static double lower_bound() { return 0.0; }
142 : };
143 :
144 0 : struct MinimumDensity {
145 0 : using type = double;
146 0 : static constexpr Options::String help = {
147 : "Minimum valid density rho_min,"
148 : " for this parametrization"};
149 0 : static double lower_bound() { return 0.0; }
150 : };
151 0 : struct MaximumDensity {
152 0 : using type = double;
153 0 : static constexpr Options::String help = {"Maximum density for this EoS"};
154 0 : static double lower_bound() { return 0.0; }
155 : };
156 :
157 0 : struct PolynomialCoefficients {
158 0 : using type = std::vector<double>;
159 0 : static constexpr Options::String help = {"Polynomial coefficients a_i"};
160 : };
161 :
162 0 : struct TrigScaling {
163 0 : using type = double;
164 0 : static constexpr Options::String help = {
165 : "Fundamental wavenumber of trig "
166 : "functions, k"};
167 0 : static double lower_bound() { return 0.0; }
168 : };
169 :
170 0 : struct SinCoefficients {
171 0 : using type = std::vector<double>;
172 0 : static constexpr Options::String help = {"Sine coefficients b_j"};
173 : };
174 0 : struct CosCoefficients {
175 0 : using type = std::vector<double>;
176 0 : static constexpr Options::String help = {"Cosine coefficients c_j"};
177 : };
178 0 : struct StitchedLowDensityEoS {
179 0 : using type = LowDensityEoS;
180 0 : static std::string name() {
181 : return pretty_type::short_name<LowDensityEoS>();
182 : }
183 0 : static constexpr Options::String help = {
184 : "Low density EoS stitched at the MinimumDensity"};
185 : };
186 :
187 0 : struct TransitionDeltaEpsilon {
188 0 : using type = double;
189 0 : static constexpr Options::String help = {
190 : "the change in internal energy across the low-"
191 : "to-high-density transition, generically 0.0"};
192 0 : static double lower_bound() { return 0.0; }
193 : };
194 :
195 0 : static constexpr Options::String help = {
196 : "An EoS with a parametrized value h(log(rho/rho_0)) with h the specific "
197 : "enthalpy and rho the baryon rest mass density. The enthalpy is "
198 : "expanded as a sum of polynomial terms and trigonometric corrections. "
199 : "let x = log(rho/rho_0) in"
200 : "h(x) = \\sum_i a_ix^i + \\sum_j b_jsin(k * j * x) + c_jcos(k * j * x) "
201 : "Note that rho(x)(1+epsilon(x)) = int_0^x e^x' h((x') dx' can be "
202 : "computed "
203 : "analytically, and therefore so can "
204 : "P(x) = rho(x) * (h(x) - (1 + epsilon(x))) "};
205 :
206 0 : using options =
207 : tmpl::list<ReferenceDensity, MaximumDensity, MinimumDensity, TrigScaling,
208 : PolynomialCoefficients, SinCoefficients, CosCoefficients,
209 : StitchedLowDensityEoS, TransitionDeltaEpsilon>;
210 :
211 0 : Enthalpy() = default;
212 0 : Enthalpy(const Enthalpy&) = default;
213 0 : Enthalpy& operator=(const Enthalpy&) = default;
214 0 : Enthalpy(Enthalpy&&) = default;
215 0 : Enthalpy& operator=(Enthalpy&&) = default;
216 0 : ~Enthalpy() override = default;
217 :
218 0 : Enthalpy(double reference_density, double max_density, double min_density,
219 : double trig_scale,
220 : const std::vector<double>& polynomial_coefficients,
221 : const std::vector<double>& sin_coefficients,
222 : const std::vector<double>& cos_coefficients,
223 : const LowDensityEoS& low_density_eos,
224 : const double transition_delta_epsilon);
225 :
226 0 : std::unique_ptr<EquationOfState<true, 1>> get_clone() const override;
227 :
228 0 : std::unique_ptr<EquationOfState<true, 3>> promote_to_3d_eos() const override;
229 :
230 0 : std::unique_ptr<EquationOfState<true, 2>> promote_to_2d_eos() const override;
231 :
232 0 : bool is_equal(const EquationOfState<true, 1>& rhs) const override;
233 :
234 0 : bool operator==(const Enthalpy<LowDensityEoS>& rhs) const;
235 :
236 0 : bool operator!=(const Enthalpy<LowDensityEoS>& rhs) const;
237 :
238 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(Enthalpy, 1)
239 :
240 0 : WRAPPED_PUPable_decl_base_template( // NOLINT
241 : SINGLE_ARG(EquationOfState<true, 1>), Enthalpy);
242 :
243 : /// The lower bound of the rest mass density that is valid for this EOS
244 1 : double rest_mass_density_lower_bound() const override { return 0.0; }
245 :
246 : /// The upper bound of the rest mass density that is valid for this EOS
247 1 : double rest_mass_density_upper_bound() const override {
248 : return std::numeric_limits<double>::max();
249 : }
250 :
251 : /// The lower bound of the specific enthalpy that is valid for this EOS
252 1 : double specific_enthalpy_lower_bound() const override { return 1.0; }
253 :
254 : /// The lower bound of the specific internal energy that is valid for this EOS
255 1 : double specific_internal_energy_lower_bound() const override { return 0.0; }
256 :
257 : /// The upper bound of the specific internal energy that is valid for this EOS
258 1 : double specific_internal_energy_upper_bound() const override {
259 : return std::numeric_limits<double>::max();
260 : }
261 :
262 : /// The vacuum baryon mass for this EoS
263 1 : double baryon_mass() const override { return low_density_eos_.baryon_mass(); }
264 :
265 : private:
266 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(1)
267 :
268 : SPECTRE_ALWAYS_INLINE
269 0 : bool in_low_density_domain(const double rest_mass_density) const {
270 : return rest_mass_density < minimum_density_;
271 : }
272 :
273 0 : double x_from_density(const double rest_mass_density) const;
274 0 : double density_from_x(const double x) const;
275 0 : double energy_density_from_log_density(const double x,
276 : const double rest_mass_density) const;
277 0 : static double evaluate_coefficients(
278 : const Enthalpy::Coefficients& coefficients, const double x,
279 : const double exponential_prefactor =
280 : std::numeric_limits<double>::signaling_NaN());
281 0 : static Enthalpy::Coefficients compute_pressure_coefficients(
282 : const typename Enthalpy::Coefficients& enthalpy,
283 : const typename Enthalpy::Coefficients& energy_density);
284 :
285 0 : double chi_from_density(const double rest_mass_density) const;
286 0 : double specific_internal_energy_from_density(
287 : const double rest_mass_density) const;
288 0 : double specific_enthalpy_from_density(const double rest_mass_density) const;
289 0 : double pressure_from_density(const double rest_mass_density) const;
290 0 : double pressure_from_log_density(const double x,
291 : const double rest_mass_density) const;
292 0 : double rest_mass_density_from_enthalpy(const double specific_enthalpy) const;
293 :
294 0 : double reference_density_ = std::numeric_limits<double>::signaling_NaN();
295 0 : double minimum_density_ = std::numeric_limits<double>::signaling_NaN();
296 0 : double maximum_density_ = std::numeric_limits<double>::signaling_NaN();
297 0 : double minimum_enthalpy_ = std::numeric_limits<double>::signaling_NaN();
298 :
299 0 : LowDensityEoS low_density_eos_;
300 0 : Coefficients coefficients_;
301 0 : Coefficients exponential_integral_coefficients_;
302 0 : Coefficients derivative_coefficients_;
303 0 : Coefficients pressure_coefficients_;
304 : };
305 :
306 : } // namespace EquationsOfState
|
