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 :
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/Serialization/CharmPupable.hpp"
22 : #include "Utilities/TMPL.hpp"
23 :
24 : /// \cond
25 : class DataVector;
26 : /// \endcond
27 :
28 : namespace EquationsOfState {
29 : /*!
30 : * \ingroup EquationsOfStateGroup
31 : *
32 : * \brief Hybrid equation of state combining a barotropic EOS for cold
33 : * (zero-temperature) part with a simple thermal part
34 : *
35 : * The hybrid equation of state:
36 : *
37 : * \f[
38 : * p = p_{cold}(\rho) + \rho (\Gamma_{th}-1) (\epsilon - \epsilon_{cold}(\rho))
39 : * \f]
40 : *
41 : * where \f$p\f$ is the pressure, \f$\rho\f$ is the rest mass density,
42 : * \f$\epsilon\f$ is the specific internal energy, \f$p_{cold}\f$ and
43 : * \f$\epsilon_{cold}\f$ are the pressure and specific internal energy evaluated
44 : * using the cold EOS, and \f$\Gamma_{th}\f$ is the adiabatic index for the
45 : * thermal part.
46 : *
47 : * The temperature \f$T\f$ is defined as
48 : *
49 : * \f[
50 : * T = (\Gamma_{th} - 1) (\epsilon - \epsilon_{cold})
51 : * \f]
52 : *
53 : * This is the amount of internal energy above that of the cold EOS.
54 : *
55 : * For a hybrid EOS with the cold EOS being a polytrope we have
56 : *
57 : * \f{align}
58 : * p &= \kappa \rho^\Gamma + \rho T
59 : * =\kappa\rho^\Gamma+\rho
60 : * (\Gamma-1)\left(\epsilon-\frac{K\rho^{\Gamma-1}}{\Gamma-1}\right), \\
61 : * T&=(\Gamma-1)\left(\epsilon-\frac{K\rho^{\Gamma-1}}{\Gamma-1}\right), \\
62 : * \epsilon &= \frac{1}{(\Gamma-1)}\frac{p}{\rho}
63 : * \f}
64 : */
65 : template <typename ColdEquationOfState>
66 1 : class HybridEos
67 : : public EquationOfState<ColdEquationOfState::is_relativistic, 2> {
68 : public:
69 0 : static constexpr size_t thermodynamic_dim = 2;
70 0 : static constexpr bool is_relativistic = ColdEquationOfState::is_relativistic;
71 :
72 0 : struct ColdEos {
73 0 : using type = ColdEquationOfState;
74 0 : static constexpr Options::String help = {"Cold equation of state"};
75 0 : static std::string name() {
76 : return pretty_type::short_name<ColdEquationOfState>();
77 : }
78 : };
79 :
80 0 : struct ThermalAdiabaticIndex {
81 0 : using type = double;
82 0 : static constexpr Options::String help = {"Adiabatic index Gamma_th"};
83 : };
84 :
85 0 : struct MinTemperature {
86 0 : using type = double;
87 0 : static type lower_bound() { return 0.0; }
88 0 : static constexpr Options::String help = {
89 : "Minimum temperature. "
90 : "This value must be non-negative."};
91 : };
92 :
93 0 : static constexpr Options::String help = {
94 : "A hybrid equation of state combining a cold EOS with a simple thermal "
95 : "part. The pressure is related to the rest mass density by "
96 : " p = p_cold(rho) + rho * (Gamma_th - 1) * (epsilon - "
97 : "epsilon_cold(rho)), where p is the pressure, rho is the rest mass "
98 : "density, epsilon is the specific internal energy, p_cold and "
99 : "epsilon_cold are the pressure and specific internal energy evaluated "
100 : "using the cold EOS and Gamma_th is the adiabatic index for the thermal "
101 : "part."};
102 :
103 0 : using options = tmpl::list<ColdEos, ThermalAdiabaticIndex, MinTemperature>;
104 :
105 0 : HybridEos() = default;
106 0 : HybridEos(const HybridEos&) = default;
107 0 : HybridEos& operator=(const HybridEos&) = default;
108 0 : HybridEos(HybridEos&&) = default;
109 0 : HybridEos& operator=(HybridEos&&) = default;
110 0 : ~HybridEos() override = default;
111 :
112 0 : HybridEos(ColdEquationOfState cold_eos, double thermal_adiabatic_index,
113 : double min_temperature = 0.0);
114 :
115 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(HybridEos, 2)
116 :
117 0 : WRAPPED_PUPable_decl_base_template( // NOLINT
118 : SINGLE_ARG(EquationOfState<is_relativistic, 2>), HybridEos);
119 :
120 0 : std::unique_ptr<EquationOfState<is_relativistic, 2>> get_clone()
121 : const override;
122 :
123 0 : std::unique_ptr<EquationOfState<is_relativistic, 3>> promote_to_3d_eos()
124 : const override;
125 :
126 : /// \brief Returns `true` if the EOS is barotropic
127 1 : bool is_barotropic() const override { return false; }
128 :
129 0 : bool operator==(const HybridEos<ColdEquationOfState>& rhs) const;
130 :
131 0 : bool operator!=(const HybridEos<ColdEquationOfState>& rhs) const;
132 :
133 0 : bool is_equal(const EquationOfState<is_relativistic, 2>& rhs) const override;
134 :
135 0 : static std::string name() {
136 : return "HybridEos(" + pretty_type::name<ColdEquationOfState>() + ")";
137 : }
138 :
139 : /// The lower bound of the rest mass density that is valid for this EOS
140 1 : double rest_mass_density_lower_bound() const override {
141 : return cold_eos_.rest_mass_density_lower_bound();
142 : }
143 :
144 : /// The upper bound of the rest mass density that is valid for this EOS
145 1 : double rest_mass_density_upper_bound() const override {
146 : return cold_eos_.rest_mass_density_upper_bound();
147 : }
148 :
149 : /// The lower bound of the specific internal energy that is valid for this EOS
150 : /// at the given rest mass density \f$\rho\f$
151 1 : double specific_internal_energy_lower_bound(
152 : const double rest_mass_density) const override {
153 : return get(cold_eos_.specific_internal_energy_from_density(
154 : Scalar<double>{rest_mass_density})) +
155 : (min_temperature_) / (thermal_adiabatic_index_ - 1.0);
156 : }
157 :
158 : /// The upper bound of the specific internal energy that is valid for this EOS
159 : /// at the given rest mass density \f$\rho\f$
160 1 : double specific_internal_energy_upper_bound(
161 : const double /*rest_mass_density*/) const override {
162 : return std::numeric_limits<double>::max();
163 : }
164 :
165 : /// The lower bound of the specific enthalpy that is valid for this EOS
166 1 : double specific_enthalpy_lower_bound() const override {
167 : return cold_eos_.specific_enthalpy_lower_bound() +
168 : (thermal_adiabatic_index_ * min_temperature_) /
169 : (thermal_adiabatic_index_ - 1.0);
170 : }
171 :
172 : /// The lower bound of the temperature that is valid for this EOS.
173 : /// Non-zero lower bound could be set to impose floor on the specific
174 : /// internal energy.
175 1 : double temperature_lower_bound() const override { return min_temperature_; }
176 :
177 : /// The vacuum baryon mass for this EoS
178 1 : double baryon_mass() const override { return cold_eos_.baryon_mass(); }
179 :
180 : private:
181 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(2)
182 :
183 0 : double specific_entropy_from_density_and_thermal_energy(
184 : double density, double energy) const;
185 :
186 0 : ColdEquationOfState cold_eos_;
187 0 : double thermal_adiabatic_index_ =
188 : std::numeric_limits<double>::signaling_NaN();
189 0 : double min_temperature_ = std::numeric_limits<double>::signaling_NaN();
190 : };
191 :
192 : /// \cond
193 : template <typename ColdEquationOfState>
194 : PUP::able::PUP_ID EquationsOfState::HybridEos<ColdEquationOfState>::my_PUP_ID =
195 : 0;
196 : /// \endcond
197 : } // namespace EquationsOfState
|
