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 :
16 : #include "DataStructures/Tensor/TypeAliases.hpp"
17 : #include "Options/String.hpp"
18 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
19 : #include "PointwiseFunctions/Hydro/Units.hpp"
20 : #include "Utilities/Serialization/CharmPupable.hpp"
21 : #include "Utilities/TMPL.hpp"
22 :
23 : /// \cond
24 : class DataVector;
25 : /// \endcond
26 :
27 : namespace EquationsOfState {
28 : /*!
29 : * \ingroup EquationsOfStateGroup
30 : * \brief Equation of state for a piecewise polytropic fluid
31 : *
32 : * A piecewise polytropic equation of state \f$p=K_i\rho^{\Gamma_i}\f$ where
33 : * \f$K_i\f$ is the polytropic constant and \f$\Gamma_i\f$ is the polytropic
34 : * exponent. Here the subscript \f$i\f$ indicates two pairs of constants and
35 : * exponents which characterize 'the stiffness' of the matter at low and high
36 : * densities. For a given density, the polytropic exponent is related to the
37 : * polytropic index \f$N_p\f$ by \f$N_p=1/(\Gamma-1)\f$. For posterity,
38 : * this two piece polytrope has been used in toy models of CCSNe (e.g.,
39 : * \cite Dimmelmeier2001 ) and could be extended to a general "M" number of
40 : * parts for simplified equations of state for neutron stars (e.g.,
41 : * \cite OBoyle2020 ). For a reference to a general piecewise polytrope, see
42 : * Section 2.4.7 of \cite RezzollaBook.
43 : */
44 : template <bool IsRelativistic>
45 1 : class PiecewisePolytropicFluid : public EquationOfState<IsRelativistic, 1> {
46 : public:
47 0 : static constexpr size_t thermodynamic_dim = 1;
48 0 : static constexpr bool is_relativistic = IsRelativistic;
49 :
50 : /// The density demarcating the high and low density descriptions of the
51 : /// fluid.
52 1 : struct PiecewisePolytropicTransitionDensity {
53 0 : using type = double;
54 0 : static constexpr Options::String help = {
55 : "Density below (above) which, the matter is described by a low (high) "
56 : "density polytropic fluid."};
57 0 : static double lower_bound() { return 0.0; }
58 : };
59 :
60 : /// The constant \f$K\f$ scaling the low density material
61 : /// \f$p=K\rho^{\Gamma}\f$.
62 : ///
63 : /// Note, by enforcing pressure continuity at the transition density
64 : /// \f$\bar{\rho}\f$, the high density constant \f$K_{high}\f$ is given
65 : /// as \f$K_{high} = K_{low} (\bar{\rho})^{\Gamma_{low} - \Gamma_{high}}\f$.
66 1 : struct PolytropicConstantLow {
67 0 : using type = double;
68 0 : static constexpr Options::String help = {
69 : "Polytropic constant K for lower"
70 : " density material"};
71 0 : static double lower_bound() { return 0.0; }
72 : };
73 :
74 : /// The exponent \f$\Gamma\f$, scaling the low density material
75 : /// \f$p=K\rho^{\Gamma}\f$.
76 1 : struct PolytropicExponentLow {
77 0 : using type = double;
78 0 : static constexpr Options::String help = {
79 : "Polytropic exponent for lower"
80 : " density material."};
81 0 : static double lower_bound() { return 1.0; }
82 : };
83 :
84 : /// The exponent \f$\Gamma\f$, scaling the high density material
85 : /// \f$p=K\rho^{\Gamma}\f$.
86 1 : struct PolytropicExponentHigh {
87 0 : using type = double;
88 0 : static constexpr Options::String help = {
89 : "Polytropic exponent for higher"
90 : " density material."};
91 0 : static double lower_bound() { return 1.0; }
92 : };
93 :
94 0 : static constexpr Options::String help = {
95 : "A piecewise polytropic fluid equation of state.\n"
96 : "The pressure is related to the rest mass density by p = K_i rho ^ "
97 : "Gamma_i, "
98 : "where p is the pressure, rho is the rest mass density, K_i is the "
99 : "polytropic constant either describing the low or high density material, "
100 : "and Gamma_i is the polytropic exponent for the low or high density "
101 : "material. The polytropic index N_i is defined as Gamma_i = 1 + 1 / N_i."
102 : " The subscript 'i' refers to different pairs of Gamma and K that can"
103 : " describe either low or high density material."};
104 :
105 0 : using options =
106 : tmpl::list<PiecewisePolytropicTransitionDensity, PolytropicConstantLow,
107 : PolytropicExponentLow, PolytropicExponentHigh>;
108 :
109 0 : PiecewisePolytropicFluid() = default;
110 0 : PiecewisePolytropicFluid(const PiecewisePolytropicFluid&) = default;
111 0 : PiecewisePolytropicFluid& operator=(const PiecewisePolytropicFluid&) =
112 : default;
113 0 : PiecewisePolytropicFluid(PiecewisePolytropicFluid&&) = default;
114 0 : PiecewisePolytropicFluid& operator=(PiecewisePolytropicFluid&&) = default;
115 0 : ~PiecewisePolytropicFluid() override = default;
116 :
117 0 : PiecewisePolytropicFluid(double transition_density,
118 : double polytropic_constant_lo,
119 : double polytropic_exponent_lo,
120 : double polytropic_exponent_hi);
121 :
122 0 : std::unique_ptr<EquationOfState<IsRelativistic, 1>> get_clone()
123 : const override;
124 :
125 1 : std::unique_ptr<EquationOfState<IsRelativistic, 3>> promote_to_3d_eos()
126 : const override;
127 :
128 1 : std::unique_ptr<EquationOfState<IsRelativistic, 2>> promote_to_2d_eos()
129 : const override;
130 :
131 0 : bool is_equal(const EquationOfState<IsRelativistic, 1>& rhs) const override;
132 :
133 0 : bool operator==(const PiecewisePolytropicFluid<IsRelativistic>& rhs) const;
134 :
135 0 : bool operator!=(const PiecewisePolytropicFluid<IsRelativistic>& rhs) const;
136 :
137 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(PiecewisePolytropicFluid, 1)
138 :
139 0 : WRAPPED_PUPable_decl_base_template( // NOLINT
140 : SINGLE_ARG(EquationOfState<IsRelativistic, 1>), PiecewisePolytropicFluid);
141 :
142 : /// The lower bound of the rest mass density that is valid for this EOS
143 1 : double rest_mass_density_lower_bound() const override { return 0.0; }
144 :
145 : /// The upper bound of the rest mass density that is valid for this EOS
146 1 : double rest_mass_density_upper_bound() const override;
147 :
148 : /// The lower bound of the specific enthalpy that is valid for this EOS
149 1 : double specific_enthalpy_lower_bound() const override {
150 : return IsRelativistic ? 1.0 : 0.0;
151 : }
152 :
153 : /// The lower bound of the specific internal energy that is valid for this EOS
154 1 : double specific_internal_energy_lower_bound() const override { return 0.0; }
155 :
156 : /// The upper bound of the specific internal energy that is valid for this EOS
157 1 : double specific_internal_energy_upper_bound() const override;
158 :
159 : /// The vacuum baryon mass for this EoS
160 1 : double baryon_mass() const override {
161 : return hydro::units::geometric::default_baryon_mass;
162 : }
163 :
164 : private:
165 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(1)
166 :
167 0 : double transition_density_ = std::numeric_limits<double>::signaling_NaN();
168 0 : double transition_pressure_ = std::numeric_limits<double>::signaling_NaN();
169 0 : double transition_spec_eint_ = std::numeric_limits<double>::signaling_NaN();
170 0 : double polytropic_constant_lo_ = std::numeric_limits<double>::signaling_NaN();
171 0 : double polytropic_exponent_lo_ = std::numeric_limits<double>::signaling_NaN();
172 0 : double polytropic_constant_hi_ = std::numeric_limits<double>::signaling_NaN();
173 0 : double polytropic_exponent_hi_ = std::numeric_limits<double>::signaling_NaN();
174 : };
175 :
176 : /// \cond
177 : template <bool IsRelativistic>
178 : PUP::able::PUP_ID
179 : EquationsOfState::PiecewisePolytropicFluid<IsRelativistic>::my_PUP_ID = 0;
180 : /// \endcond
181 : } // namespace EquationsOfState
|
