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 polytropic fluid
31 : *
32 : * A polytropic equation of state \f$p=K\rho^{\Gamma}\f$ where \f$K\f$ is the
33 : * polytropic constant and \f$\Gamma\f$ is the polytropic exponent. The
34 : * polytropic exponent is related to the polytropic index \f$N_p\f$ by
35 : * \f$N_p=1/(\Gamma-1)\f$.
36 : *
37 : * We also have
38 : *
39 : * \f{align}{
40 : * \epsilon&=\frac{K\rho^{\Gamma-1}}{\Gamma-1}\\
41 : * h&=1+\epsilon+\frac{p}{\rho}=1+\frac{K\Gamma}{\Gamma-1}\rho^{\Gamma-1} \\
42 : * T&=0 \\
43 : * c_s^2&=\frac{\Gamma p}{\rho h}
44 : * =\frac{\Gamma(\Gamma-1)p}{\rho(\Gamma-1)+\Gamma p}
45 : * =\left(\frac{1}{\Gamma K\rho^{\Gamma-1}}+\frac{1}{\Gamma-1}\right)^{-1}
46 : * \f}
47 : */
48 : template <bool IsRelativistic>
49 1 : class PolytropicFluid : public EquationOfState<IsRelativistic, 1> {
50 : public:
51 0 : static constexpr size_t thermodynamic_dim = 1;
52 0 : static constexpr bool is_relativistic = IsRelativistic;
53 :
54 0 : struct PolytropicConstant {
55 0 : using type = double;
56 0 : static constexpr Options::String help = {"Polytropic constant K"};
57 0 : static double lower_bound() { return 0.0; }
58 : };
59 :
60 0 : struct PolytropicExponent {
61 0 : using type = double;
62 0 : static constexpr Options::String help = {"Polytropic exponent Gamma"};
63 0 : static double lower_bound() { return 1.0; }
64 : };
65 :
66 0 : static constexpr Options::String help = {
67 : "A polytropic fluid equation of state.\n"
68 : "The pressure is related to the rest mass density by p = K rho ^ Gamma, "
69 : "where p is the pressure, rho is the rest mass density, K is the "
70 : "polytropic constant, and Gamma is the polytropic exponent. The "
71 : "polytropic index N is defined as Gamma = 1 + 1 / N."};
72 :
73 0 : using options = tmpl::list<PolytropicConstant, PolytropicExponent>;
74 :
75 0 : PolytropicFluid() = default;
76 0 : PolytropicFluid(const PolytropicFluid&) = default;
77 0 : PolytropicFluid& operator=(const PolytropicFluid&) = default;
78 0 : PolytropicFluid(PolytropicFluid&&) = default;
79 0 : PolytropicFluid& operator=(PolytropicFluid&&) = default;
80 0 : ~PolytropicFluid() override = default;
81 :
82 0 : PolytropicFluid(double polytropic_constant, double polytropic_exponent);
83 :
84 0 : std::unique_ptr<EquationOfState<IsRelativistic, 1>> get_clone()
85 : const override;
86 :
87 1 : std::unique_ptr<EquationOfState<IsRelativistic, 3>> promote_to_3d_eos()
88 : const override;
89 :
90 1 : std::unique_ptr<EquationOfState<IsRelativistic, 2>> promote_to_2d_eos()
91 : const override;
92 :
93 0 : bool is_equal(const EquationOfState<IsRelativistic, 1>& rhs) const override;
94 :
95 0 : bool operator==(const PolytropicFluid<IsRelativistic>& rhs) const;
96 :
97 0 : bool operator!=(const PolytropicFluid<IsRelativistic>& rhs) const;
98 :
99 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(PolytropicFluid, 1)
100 :
101 0 : WRAPPED_PUPable_decl_base_template( // NOLINT
102 : SINGLE_ARG(EquationOfState<IsRelativistic, 1>), PolytropicFluid);
103 :
104 : /// The lower bound of the rest mass density that is valid for this EOS
105 1 : double rest_mass_density_lower_bound() const override { return 0.0; }
106 :
107 : /// The upper bound of the rest mass density that is valid for this EOS
108 1 : double rest_mass_density_upper_bound() const override;
109 :
110 : /// The lower bound of the specific enthalpy that is valid for this EOS
111 1 : double specific_enthalpy_lower_bound() const override {
112 : return IsRelativistic ? 1.0 : 0.0;
113 : }
114 :
115 : /// The lower bound of the specific internal energy that is valid for this EOS
116 1 : double specific_internal_energy_lower_bound() const override { return 0.0; }
117 :
118 : /// The upper bound of the specific internal energy that is valid for this EOS
119 1 : double specific_internal_energy_upper_bound() const override;
120 :
121 : /// The vacuum baryon mass for this EoS
122 1 : double baryon_mass() const override {
123 : return hydro::units::geometric::default_baryon_mass;
124 : }
125 :
126 : private:
127 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(1)
128 :
129 0 : double polytropic_constant_ = std::numeric_limits<double>::signaling_NaN();
130 0 : double polytropic_exponent_ = std::numeric_limits<double>::signaling_NaN();
131 : };
132 :
133 : /// \cond
134 : template <bool IsRelativistic>
135 : PUP::able::PUP_ID EquationsOfState::PolytropicFluid<IsRelativistic>::my_PUP_ID =
136 : 0;
137 : /// \endcond
138 : } // namespace EquationsOfState
|
