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/Tensor.hpp"
17 : #include "DataStructures/Tensor/TypeAliases.hpp"
18 : #include "IO/H5/EosTable.hpp"
19 : #include "IO/H5/File.hpp"
20 : #include "NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp"
21 : #include "Options/String.hpp"
22 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
23 : #include "PointwiseFunctions/Hydro/Units.hpp"
24 : #include "Utilities/Serialization/CharmPupable.hpp"
25 : #include "Utilities/TMPL.hpp"
26 :
27 : /// \cond
28 : class DataVector;
29 : /// \endcond
30 :
31 : namespace EquationsOfState {
32 : /*!
33 : * \ingroup EquationsOfStateGroup
34 : * \brief Nuclear matter equation of state in tabulated form.
35 : *
36 : * The equation of state takes the form
37 : *
38 : * \f[
39 : * p = p (T, rho, Y_e)
40 : * \f]
41 : *
42 : * where \f$\rho\f$ is the rest mass density, \f$T\f$ is the
43 : * temperature, and \f$Y_e\f$ is the electron fraction.
44 : * The temperature is given in units of MeV.
45 : */
46 : template <bool IsRelativistic>
47 1 : class Tabulated3D : public EquationOfState<IsRelativistic, 3> {
48 : public:
49 0 : static constexpr size_t thermodynamic_dim = 3;
50 0 : static constexpr bool is_relativistic = IsRelativistic;
51 :
52 0 : static constexpr Options::String help = {
53 : "A tabulated three-dimensional equation of state.\n"
54 : "The energy density, pressure and sound speed "
55 : "are tabulated as a function of density, electron_fraction and "
56 : "temperature."};
57 :
58 0 : struct TableFilename {
59 0 : using type = std::string;
60 0 : static constexpr Options::String help{"File name of the EOS table"};
61 : };
62 :
63 0 : struct TableSubFilename {
64 0 : using type = std::string;
65 0 : static constexpr Options::String help{
66 : "Subfile name of the EOS table, e.g., 'dd2'."};
67 : };
68 :
69 0 : using options = tmpl::list<TableFilename, TableSubFilename>;
70 :
71 : /// Fields stored in the table
72 0 : enum : size_t {
73 : Epsilon = 0,
74 : Pressure,
75 : CsSquared,
76 : DeltaMu,
77 : Kappa,
78 : Zeta,
79 : SpecificEntropy,
80 : NumberOfVars
81 : };
82 :
83 0 : Tabulated3D() = default;
84 0 : Tabulated3D(const Tabulated3D&) = default;
85 0 : Tabulated3D& operator=(const Tabulated3D&) = default;
86 0 : Tabulated3D(Tabulated3D&&) = default;
87 0 : Tabulated3D& operator=(Tabulated3D&&) = default;
88 0 : ~Tabulated3D() override = default;
89 :
90 0 : explicit Tabulated3D(const std::string& filename,
91 : const std::string& subfilename);
92 :
93 0 : explicit Tabulated3D(std::vector<double> electron_fraction,
94 : std::vector<double> log_density,
95 : std::vector<double> log_temperature,
96 : std::vector<double> table_data, double energy_shift,
97 : double enthalpy_minimum);
98 :
99 0 : explicit Tabulated3D(const h5::EosTable& spectre_eos);
100 :
101 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(Tabulated3D, 3)
102 :
103 : template <class DataType>
104 0 : void convert_to_table_quantities(
105 : const gsl::not_null<Scalar<DataType>*> converted_electron_fraction,
106 : const gsl::not_null<Scalar<DataType>*> log_rest_mass_density,
107 : const gsl::not_null<Scalar<DataType>*> log_temperature,
108 : const Scalar<DataType>& electron_fraction,
109 : const Scalar<DataType>& rest_mass_density,
110 : const Scalar<DataType>& temperature) const {
111 : get(*converted_electron_fraction) = get(electron_fraction);
112 : get(*log_rest_mass_density) = get(rest_mass_density);
113 : get(*log_temperature) = get(temperature);
114 :
115 : // Enforce physicality of input
116 : // We reuse the same variables here, these are not log yet.
117 : enforce_physicality(*converted_electron_fraction, *log_rest_mass_density,
118 : *log_temperature);
119 :
120 : // Table uses log T and log rho
121 : get(*log_rest_mass_density) = log(get(*log_rest_mass_density));
122 : get(*log_temperature) = log(get(*log_temperature));
123 : }
124 :
125 0 : std::unique_ptr<EquationOfState<IsRelativistic, 3>> get_clone()
126 : const override;
127 :
128 0 : void initialize(std::vector<double> electron_fraction,
129 : std::vector<double> log_density,
130 : std::vector<double> log_temperature,
131 : std::vector<double> table_data, double energy_shift,
132 : double enthalpy_minimum);
133 :
134 :
135 0 : void initialize(const h5::EosTable& spectre_eos);
136 :
137 0 : bool is_equal(const EquationOfState<IsRelativistic, 3>& rhs) const override;
138 :
139 : /// \brief Returns `true` if the EOS is barotropic
140 1 : bool is_barotropic() const override { return false; }
141 :
142 : /// \brief Returns `true` if the EOS is in beta-equilibrium
143 1 : bool is_equilibrium() const override { return false; }
144 :
145 0 : bool operator==(const Tabulated3D<IsRelativistic>& rhs) const;
146 :
147 0 : bool operator!=(const Tabulated3D<IsRelativistic>& rhs) const;
148 :
149 : template <class DataType>
150 0 : Scalar<DataType> equilibrium_electron_fraction_from_density_temperature_impl(
151 : const Scalar<DataType>& rest_mass_density,
152 : const Scalar<DataType>& temperature) const;
153 :
154 : /// @{
155 : /*!
156 : * Computes the electron fraction in beta-equilibrium \f$Y_e^{\rm eq}\f$ from
157 : * the rest mass density \f$\rho\f$ and the temperature \f$T\f$.
158 : */
159 1 : Scalar<double> equilibrium_electron_fraction_from_density_temperature(
160 : const Scalar<double>& rest_mass_density,
161 : const Scalar<double>& temperature) const override {
162 : return equilibrium_electron_fraction_from_density_temperature_impl<double>(
163 : rest_mass_density, temperature);
164 : }
165 :
166 1 : Scalar<DataVector> equilibrium_electron_fraction_from_density_temperature(
167 : const Scalar<DataVector>& rest_mass_density,
168 : const Scalar<DataVector>& temperature) const override {
169 : return equilibrium_electron_fraction_from_density_temperature_impl<
170 : DataVector>(rest_mass_density, temperature);
171 : }
172 : /// @}
173 : //
174 :
175 : template <typename DataType>
176 0 : void enforce_physicality(Scalar<DataType>& electron_fraction,
177 : Scalar<DataType>& density,
178 : Scalar<DataType>& temperature) const;
179 :
180 0 : WRAPPED_PUPable_decl_base_template( // NOLINT
181 : SINGLE_ARG(EquationOfState<IsRelativistic, 3>), Tabulated3D);
182 :
183 : /// The lower bound of the electron fraction that is valid for this EOS
184 1 : double electron_fraction_lower_bound() const override {
185 : return table_electron_fraction_.front();
186 : }
187 :
188 : /// The upper bound of the electron fraction that is valid for this EOS
189 1 : double electron_fraction_upper_bound() const override {
190 : return table_electron_fraction_.back();
191 : }
192 :
193 : /// The lower bound of the rest mass density that is valid for this EOS
194 1 : double rest_mass_density_lower_bound() const override {
195 : return std::exp((table_log_density_.front()));
196 : }
197 :
198 : /// The upper bound of the rest mass density that is valid for this EOS
199 1 : double rest_mass_density_upper_bound() const override {
200 : return std::exp((table_log_density_.back()));
201 : }
202 :
203 : /// The lower bound of the temperature that is valid for this EOS
204 1 : double temperature_lower_bound() const override {
205 : return std::exp((table_log_temperature_.front()));
206 : }
207 :
208 : /// The upper bound of the temperature that is valid for this EOS
209 1 : double temperature_upper_bound() const override {
210 : return std::exp((table_log_temperature_.back()));
211 : }
212 :
213 : /// The lower bound of the specific internal energy that is valid for this EOS
214 : /// at the given rest mass density \f$\rho\f$ and electron fraction \f$Y_e\f$
215 1 : double specific_internal_energy_lower_bound(
216 : const double rest_mass_density,
217 : const double electron_fraction) const override;
218 :
219 : /// The upper bound of the specific internal energy that is valid for this EOS
220 : /// at the given rest mass density \f$\rho\f$
221 1 : double specific_internal_energy_upper_bound(
222 : const double rest_mass_density,
223 : const double electron_fraction) const override;
224 :
225 : /// The lower bound of the specific enthalpy that is valid for this EOS
226 1 : double specific_enthalpy_lower_bound() const override {
227 : return enthalpy_minimum_;
228 : }
229 :
230 : /// The baryon mass for this EoS
231 1 : double baryon_mass() const override {
232 : return hydro::units::geometric::neutron_mass;
233 : }
234 :
235 : private:
236 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(3)
237 :
238 0 : void initialize_interpolator();
239 :
240 : /// Energy shift used to account for negative specific internal energies,
241 : /// which are only stored logarithmically
242 1 : double energy_shift_ = 0.;
243 :
244 : /// Enthalpy minium across the table
245 1 : double enthalpy_minimum_ = 1.;
246 :
247 : /// Main interpolator for the EoS.
248 : /// The ordering is \f$(\log T. \log \rho, Y_e)\f$.
249 : /// Assumed to be sorted in ascending order.
250 1 : intrp::UniformMultiLinearSpanInterpolation<3, NumberOfVars> interpolator_{};
251 : /// Electron fraction
252 1 : std::vector<double> table_electron_fraction_{};
253 : /// Logarithmic rest-mass denisty
254 1 : std::vector<double> table_log_density_{};
255 : /// Logarithmic temperature
256 1 : std::vector<double> table_log_temperature_{};
257 : /// Tabulate data. Entries are stated in the enum
258 1 : std::vector<double> table_data_{};
259 :
260 : /// Tolerance on upper bound for root finding
261 1 : static constexpr double upper_bound_tolerance_ = 0.9999;
262 : };
263 :
264 : /// \cond
265 : template <bool IsRelativistic>
266 : PUP::able::PUP_ID EquationsOfState::Tabulated3D<IsRelativistic>::my_PUP_ID = 0;
267 : /// \endcond
268 :
269 : } // namespace EquationsOfState
|
