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/sub.hpp>
7 : #include <boost/preprocessor/list/for_each.hpp>
8 : #include <boost/preprocessor/punctuation/comma_if.hpp>
9 : #include <boost/preprocessor/repetition/repeat.hpp>
10 : #include <boost/preprocessor/tuple/enum.hpp>
11 : #include <boost/preprocessor/tuple/to_list.hpp>
12 :
13 : #include "DataStructures/DataVector.hpp"
14 : #include "DataStructures/Tensor/Tensor.hpp"
15 : #include "PointwiseFunctions/Hydro/Units.hpp"
16 : #include "Utilities/CallWithDynamicType.hpp"
17 : #include "Utilities/Serialization/CharmPupable.hpp"
18 : #include "Utilities/TMPL.hpp"
19 : #include "Utilities/TypeTraits.hpp"
20 :
21 : /// \cond
22 : namespace EquationsOfState {
23 : template <typename ColdEos>
24 : class Barotropic2D;
25 : template <typename ColdEquilEos>
26 : class Barotropic3D;
27 : template <bool IsRelativistic>
28 : class DarkEnergyFluid;
29 : template <typename EquilEos>
30 : class Equilibrium3D;
31 : template <typename ColdEquationOfState>
32 : class HybridEos;
33 : template <bool IsRelativistic>
34 : class IdealFluid;
35 : template <bool IsRelativistic>
36 : class PolytropicFluid;
37 : template <bool IsRelativistic>
38 : class PiecewisePolytropicFluid;
39 : class Spectral;
40 : template <typename LowDensityEoS>
41 : class Enthalpy;
42 : template <bool IsRelativistic>
43 : class Tabulated3D;
44 : } // namespace EquationsOfState
45 : /// \endcond
46 :
47 : /// Contains all equations of state, including base class
48 : namespace EquationsOfState {
49 :
50 : namespace detail {
51 : template <bool IsRelativistic, size_t ThermodynamicDim>
52 : struct DerivedClasses {};
53 :
54 : template <>
55 : struct DerivedClasses<true, 1> {
56 : using type = tmpl::list<
57 : Enthalpy<Enthalpy<Enthalpy<Spectral>>>, Enthalpy<Enthalpy<Spectral>>,
58 : Enthalpy<Spectral>, Enthalpy<PolytropicFluid<true>>,
59 : PiecewisePolytropicFluid<true>, PolytropicFluid<true>, Spectral>;
60 : };
61 :
62 : template <>
63 : struct DerivedClasses<false, 1> {
64 : using type =
65 : tmpl::list<PiecewisePolytropicFluid<false>, PolytropicFluid<false>>;
66 : };
67 :
68 : template <>
69 : struct DerivedClasses<true, 2> {
70 : using type =
71 : tmpl::list<Barotropic2D<PolytropicFluid<true>>, Barotropic2D<Spectral>,
72 : Barotropic2D<Enthalpy<Spectral>>,
73 : Barotropic2D<PiecewisePolytropicFluid<true>>,
74 : Barotropic2D<Enthalpy<Enthalpy<Spectral>>>,
75 : Barotropic2D<Enthalpy<Enthalpy<Enthalpy<Spectral>>>>,
76 : DarkEnergyFluid<true>, IdealFluid<true>,
77 : HybridEos<PolytropicFluid<true>>, HybridEos<Spectral>,
78 : HybridEos<Enthalpy<Spectral>>>;
79 : };
80 :
81 : template <>
82 : struct DerivedClasses<false, 2> {
83 : using type = tmpl::list<Barotropic2D<PolytropicFluid<false>>,
84 : Barotropic2D<PiecewisePolytropicFluid<false>>,
85 : IdealFluid<false>, HybridEos<PolytropicFluid<false>>>;
86 : };
87 :
88 : template <>
89 : struct DerivedClasses<true, 3> {
90 : using type =
91 : tmpl::list<Tabulated3D<true>, Barotropic3D<PolytropicFluid<true>>,
92 : Barotropic3D<Spectral>, Barotropic3D<Enthalpy<Spectral>>,
93 : Barotropic3D<PiecewisePolytropicFluid<true>>,
94 : Barotropic3D<Enthalpy<Enthalpy<Spectral>>>,
95 : Barotropic3D<Enthalpy<Enthalpy<Enthalpy<Spectral>>>>,
96 : Equilibrium3D<HybridEos<PolytropicFluid<true>>>,
97 : Equilibrium3D<HybridEos<Spectral>>,
98 : Equilibrium3D<HybridEos<Enthalpy<Spectral>>>,
99 : Equilibrium3D<DarkEnergyFluid<true>>,
100 : Equilibrium3D<IdealFluid<true>>>;
101 : };
102 :
103 : template <>
104 : struct DerivedClasses<false, 3> {
105 : using type = tmpl::list<Tabulated3D<false>, Equilibrium3D<IdealFluid<false>>,
106 : Barotropic3D<PiecewisePolytropicFluid<false>>,
107 : Equilibrium3D<HybridEos<PolytropicFluid<false>>>,
108 : Barotropic3D<PolytropicFluid<false>>>;
109 : };
110 :
111 : } // namespace detail
112 :
113 : /*!
114 : * \ingroup EquationsOfStateGroup
115 : * \brief Base class for equations of state depending on whether or not the
116 : * system is relativistic, and the number of independent thermodynamic variables
117 : * (`ThermodynamicDim`) needed to determine the pressure.
118 : *
119 : * The template parameter `IsRelativistic` is `true` for relativistic equations
120 : * of state and `false` for non-relativistic equations of state.
121 : */
122 : template <bool IsRelativistic, size_t ThermodynamicDim>
123 1 : class EquationOfState;
124 :
125 : template <typename T>
126 0 : struct get_eos_base_impl {
127 0 : using type = EquationsOfState::EquationOfState<T::is_relativistic,
128 : T::thermodynamic_dim>;
129 : };
130 :
131 : template <bool IsRelativistic, size_t ThermodynamicDim>
132 0 : struct get_eos_base_impl<
133 : EquationsOfState::EquationOfState<IsRelativistic, ThermodynamicDim>> {
134 0 : using type =
135 : EquationsOfState::EquationOfState<IsRelativistic, ThermodynamicDim>;
136 : };
137 :
138 : template <typename T>
139 0 : using get_eos_base = typename get_eos_base_impl<T>::type;
140 :
141 : /*!
142 : * \ingroup EquationsOfStateGroup
143 : * \brief Base class for equations of state which need one thermodynamic
144 : * variable in order to determine the pressure.
145 : *
146 : * The template parameter `IsRelativistic` is `true` for relativistic equations
147 : * of state and `false` for non-relativistic equations of state.
148 : */
149 : template <bool IsRelativistic>
150 1 : class EquationOfState<IsRelativistic, 1> : public PUP::able {
151 : public:
152 0 : static constexpr bool is_relativistic = IsRelativistic;
153 0 : static constexpr size_t thermodynamic_dim = 1;
154 0 : using creatable_classes =
155 : typename detail::DerivedClasses<IsRelativistic, 1>::type;
156 :
157 0 : EquationOfState() = default;
158 0 : EquationOfState(const EquationOfState&) = default;
159 0 : EquationOfState& operator=(const EquationOfState&) = default;
160 0 : EquationOfState(EquationOfState&&) = default;
161 0 : EquationOfState& operator=(EquationOfState&&) = default;
162 0 : ~EquationOfState() override = default;
163 :
164 0 : explicit EquationOfState(CkMigrateMessage* msg) : PUP::able(msg) {}
165 :
166 0 : WRAPPED_PUPable_abstract(EquationOfState); // NOLINT
167 :
168 0 : virtual std::unique_ptr<EquationOfState<IsRelativistic, 1>> get_clone()
169 : const = 0;
170 :
171 0 : virtual bool is_equal(
172 : const EquationOfState<IsRelativistic, 1>& rhs) const = 0;
173 :
174 : /// Create a 3D EOS from the 1D EOS
175 : virtual std::unique_ptr<EquationOfState<IsRelativistic, 3>>
176 1 : promote_to_3d_eos() const = 0;
177 :
178 : /// Create a 2D EOS from the 1D EOS
179 : virtual std::unique_ptr<EquationOfState<IsRelativistic, 2>>
180 1 : promote_to_2d_eos() const = 0;
181 :
182 : /// \brief Returns `true` if the EOS is barotropic
183 1 : bool is_barotropic() const { return true; }
184 : /// @{
185 : /*!
186 : * Computes the electron fraction in beta-equilibrium \f$Y_e^{\rm eq}\f$ from
187 : * the rest mass density \f$\rho\f$.
188 : */
189 1 : virtual Scalar<double> equilibrium_electron_fraction_from_density_temperature(
190 : const Scalar<double>& rest_mass_density,
191 : const Scalar<double>& /*temperature*/) const {
192 : return make_with_value<Scalar<double>>(rest_mass_density, 0.1);
193 : }
194 :
195 : virtual Scalar<DataVector>
196 1 : equilibrium_electron_fraction_from_density_temperature(
197 : const Scalar<DataVector>& rest_mass_density,
198 : const Scalar<DataVector>& /*temperature*/) const {
199 : return make_with_value<Scalar<DataVector>>(rest_mass_density, 0.1);
200 : }
201 : /// @}
202 :
203 : /// @{
204 : /*!
205 : * Computes the pressure \f$p\f$ from the rest mass density \f$\rho\f$.
206 : */
207 1 : virtual Scalar<double> pressure_from_density(
208 : const Scalar<double>& /*rest_mass_density*/) const = 0;
209 1 : virtual Scalar<DataVector> pressure_from_density(
210 : const Scalar<DataVector>& /*rest_mass_density*/) const = 0;
211 : /// @}
212 :
213 : /// @{
214 : /*!
215 : * Computes the rest mass density \f$\rho\f$ from the specific enthalpy
216 : * \f$h\f$.
217 : */
218 1 : virtual Scalar<double> rest_mass_density_from_enthalpy(
219 : const Scalar<double>& /*specific_enthalpy*/) const = 0;
220 1 : virtual Scalar<DataVector> rest_mass_density_from_enthalpy(
221 : const Scalar<DataVector>& /*specific_enthalpy*/) const = 0;
222 : /// @}
223 :
224 : /// @{
225 : /*!
226 : * Computes the specific internal energy \f$\epsilon\f$ from the rest mass
227 : * density \f$\rho\f$.
228 : */
229 1 : virtual Scalar<double> specific_internal_energy_from_density(
230 : const Scalar<double>& /*rest_mass_density*/) const = 0;
231 1 : virtual Scalar<DataVector> specific_internal_energy_from_density(
232 : const Scalar<DataVector>& /*rest_mass_density*/) const = 0;
233 : /// @}
234 :
235 : /// @{
236 : /*!
237 : * Computes the temperature \f$T\f$ from the rest mass
238 : * density \f$\rho\f$.
239 : */
240 1 : virtual Scalar<double> temperature_from_density(
241 : const Scalar<double>& /*rest_mass_density*/) const {
242 : return Scalar<double>{0.0};
243 : }
244 1 : virtual Scalar<DataVector> temperature_from_density(
245 : const Scalar<DataVector>& rest_mass_density) const {
246 : return make_with_value<Scalar<DataVector>>(rest_mass_density, 0.0);
247 : }
248 : /// @}
249 :
250 : /// @{
251 : /*!
252 : * Computes the temperature \f$\T\f$ from the specific internal energy
253 : * \f$\epsilon\f$.
254 : */
255 1 : virtual Scalar<double> temperature_from_specific_internal_energy(
256 : const Scalar<double>& /*specific_internal_energy*/) const {
257 : return Scalar<double>{0.0};
258 : }
259 1 : virtual Scalar<DataVector> temperature_from_specific_internal_energy(
260 : const Scalar<DataVector>& specific_internal_energy) const {
261 : return make_with_value<Scalar<DataVector>>(specific_internal_energy, 0.0);
262 : }
263 : /// @}
264 :
265 : /// @{
266 : /*!
267 : * Computes \f$\chi=\partial p / \partial \rho\f$ from \f$\rho\f$, where
268 : * \f$p\f$ is the pressure and \f$\rho\f$ is the rest mass density.
269 : */
270 1 : virtual Scalar<double> chi_from_density(
271 : const Scalar<double>& /*rest_mass_density*/) const = 0;
272 1 : virtual Scalar<DataVector> chi_from_density(
273 : const Scalar<DataVector>& /*rest_mass_density*/) const = 0;
274 : /// @}
275 :
276 : /// @{
277 : /*!
278 : * Computes \f$\kappa p/\rho^2=(p/\rho^2)\partial p / \partial \epsilon\f$
279 : * from \f$\rho\f$, where \f$p\f$ is the pressure, \f$\rho\f$ is the rest mass
280 : * density, and \f$\epsilon\f$ is the specific internal energy.
281 : *
282 : * The reason for not returning just
283 : * \f$\kappa=\partial p / \partial \epsilon\f$ is to avoid division by zero
284 : * for small values of \f$\rho\f$ when assembling the speed of sound with
285 : * some equations of state.
286 : */
287 1 : virtual Scalar<double> kappa_times_p_over_rho_squared_from_density(
288 : const Scalar<double>& /*rest_mass_density*/) const = 0;
289 1 : virtual Scalar<DataVector> kappa_times_p_over_rho_squared_from_density(
290 : const Scalar<DataVector>& /*rest_mass_density*/) const = 0;
291 :
292 : /// The lower bound of the electron fraction that is valid for this EOS
293 1 : virtual double electron_fraction_lower_bound() const { return 0.0; }
294 :
295 : /// The upper bound of the electron fraction that is valid for this EOS
296 1 : virtual double electron_fraction_upper_bound() const { return 1.0; }
297 :
298 : /// The lower bound of the rest mass density that is valid for this EOS
299 1 : virtual double rest_mass_density_lower_bound() const = 0;
300 :
301 : /// The upper bound of the rest mass density that is valid for this EOS
302 1 : virtual double rest_mass_density_upper_bound() const = 0;
303 :
304 : /// The lower bound of the temperature that is valid for this EOS
305 1 : virtual double temperature_lower_bound() const { return 0.0; };
306 :
307 : /// The upper bound of the temperature that is valid for this EOS
308 1 : virtual double temperature_upper_bound() const {
309 : return std::numeric_limits<double>::max();
310 : };
311 :
312 : /// The lower bound of the specific internal energy that is valid for this EOS
313 1 : virtual double specific_internal_energy_lower_bound() const { return 0.0; };
314 :
315 : /// The upper bound of the specific internal energy that is valid for this EOS
316 1 : virtual double specific_internal_energy_upper_bound() const {
317 : return std::numeric_limits<double>::max();
318 : };
319 :
320 : /// The lower bound of the specific enthalpy that is valid for this EOS
321 1 : virtual double specific_enthalpy_lower_bound() const = 0;
322 :
323 : /// The vacuum mass of a baryon for this EOS
324 1 : virtual double baryon_mass() const {
325 : return hydro::units::geometric::default_baryon_mass;
326 : }
327 : };
328 :
329 : /*!
330 : * \ingroup EquationsOfStateGroup
331 : * \brief Base class for equations of state which need two independent
332 : * thermodynamic variables in order to determine the pressure.
333 : *
334 : * The template parameter `IsRelativistic` is `true` for relativistic equations
335 : * of state and `false` for non-relativistic equations of state.
336 : */
337 : template <bool IsRelativistic>
338 1 : class EquationOfState<IsRelativistic, 2> : public PUP::able {
339 : public:
340 0 : static constexpr bool is_relativistic = IsRelativistic;
341 0 : static constexpr size_t thermodynamic_dim = 2;
342 0 : using creatable_classes =
343 : typename detail::DerivedClasses<IsRelativistic, 2>::type;
344 :
345 0 : EquationOfState() = default;
346 0 : EquationOfState(const EquationOfState&) = default;
347 0 : EquationOfState& operator=(const EquationOfState&) = default;
348 0 : EquationOfState(EquationOfState&&) = default;
349 0 : EquationOfState& operator=(EquationOfState&&) = default;
350 0 : ~EquationOfState() override = default;
351 :
352 0 : explicit EquationOfState(CkMigrateMessage* msg) : PUP::able(msg) {}
353 :
354 0 : WRAPPED_PUPable_abstract(EquationOfState); // NOLINT
355 :
356 0 : virtual inline std::unique_ptr<EquationOfState<IsRelativistic, 2>> get_clone()
357 : const = 0;
358 :
359 0 : virtual bool is_equal(
360 : const EquationOfState<IsRelativistic, 2>& rhs) const = 0;
361 :
362 : virtual std::unique_ptr<EquationOfState<IsRelativistic, 2>>
363 0 : promote_to_2d_eos() const {
364 : return this->get_clone();
365 : }
366 :
367 : virtual std::unique_ptr<EquationOfState<IsRelativistic, 3>>
368 0 : promote_to_3d_eos() const = 0;
369 :
370 : /// \brief Returns `true` if the EOS is barotropic
371 1 : virtual bool is_barotropic() const = 0;
372 :
373 : /// \brief Returns `true` if the EOS is in beta-equilibrium
374 1 : virtual bool is_equilibrium() const {
375 : return true;
376 : }
377 :
378 : /// @{
379 : /*!
380 : * Computes the electron fraction in beta-equilibrium \f$Y_e^{\rm eq}\f$ from
381 : * the rest mass density \f$\rho\f$ and the temperature \f$T\f$.
382 : */
383 1 : virtual Scalar<double> equilibrium_electron_fraction_from_density_temperature(
384 : const Scalar<double>& rest_mass_density,
385 : const Scalar<double>& /*temperature*/) const {
386 : return make_with_value<Scalar<double>>(rest_mass_density, 0.1);
387 : }
388 :
389 : virtual Scalar<DataVector>
390 1 : equilibrium_electron_fraction_from_density_temperature(
391 : const Scalar<DataVector>& rest_mass_density,
392 : const Scalar<DataVector>& /*temperature*/) const {
393 : return make_with_value<Scalar<DataVector>>(rest_mass_density, 0.1);
394 : }
395 : /// @}
396 :
397 : /// @{
398 : /*!
399 : * Computes the pressure \f$p\f$ from the rest mass density \f$\rho\f$ and the
400 : * specific internal energy \f$\epsilon\f$.
401 : */
402 1 : virtual Scalar<double> pressure_from_density_and_energy(
403 : const Scalar<double>& /*rest_mass_density*/,
404 : const Scalar<double>& /*specific_internal_energy*/) const = 0;
405 1 : virtual Scalar<DataVector> pressure_from_density_and_energy(
406 : const Scalar<DataVector>& /*rest_mass_density*/,
407 : const Scalar<DataVector>& /*specific_internal_energy*/) const = 0;
408 : /// @}
409 :
410 : /// @{
411 : /*!
412 : * Computes the pressure \f$p\f$ from the rest mass density \f$\rho\f$ and the
413 : * specific enthalpy \f$h\f$.
414 : */
415 1 : virtual Scalar<double> pressure_from_density_and_enthalpy(
416 : const Scalar<double>& /*rest_mass_density*/,
417 : const Scalar<double>& /*specific_enthalpy*/) const = 0;
418 1 : virtual Scalar<DataVector> pressure_from_density_and_enthalpy(
419 : const Scalar<DataVector>& /*rest_mass_density*/,
420 : const Scalar<DataVector>& /*specific_enthalpy*/) const = 0;
421 : /// @}
422 :
423 : /// @{
424 : /*!
425 : * Computes the specific internal energy \f$\epsilon\f$ from the rest mass
426 : * density \f$\rho\f$ and the pressure \f$p\f$.
427 : */
428 1 : virtual Scalar<double> specific_internal_energy_from_density_and_pressure(
429 : const Scalar<double>& /*rest_mass_density*/,
430 : const Scalar<double>& /*pressure*/) const = 0;
431 1 : virtual Scalar<DataVector> specific_internal_energy_from_density_and_pressure(
432 : const Scalar<DataVector>& /*rest_mass_density*/,
433 : const Scalar<DataVector>& /*pressure*/) const = 0;
434 : /// @}
435 :
436 : /// @{
437 : /*!
438 : * Computes the temperature \f$T\f$ from the rest mass
439 : * density \f$\rho\f$ and the specific internal energy \f$\epsilon\f$.
440 : */
441 1 : virtual Scalar<double> temperature_from_density_and_energy(
442 : const Scalar<double>& /*rest_mass_density*/,
443 : const Scalar<double>& /*specific_internal_energy*/) const = 0;
444 1 : virtual Scalar<DataVector> temperature_from_density_and_energy(
445 : const Scalar<DataVector>& /*rest_mass_density*/,
446 : const Scalar<DataVector>& /*specific_internal_energy*/) const = 0;
447 : /// @}
448 :
449 : /// @{
450 : /*!
451 : * Computes the specific internal energy \f$\epsilon\f$ from the rest mass
452 : * density \f$\rho\f$ and the temperature \f$T\f$.
453 : */
454 1 : virtual Scalar<double> specific_internal_energy_from_density_and_temperature(
455 : const Scalar<double>& /*rest_mass_density*/,
456 : const Scalar<double>& /*temperature*/) const = 0;
457 : virtual Scalar<DataVector>
458 1 : specific_internal_energy_from_density_and_temperature(
459 : const Scalar<DataVector>& /*rest_mass_density*/,
460 : const Scalar<DataVector>& /*temperature*/) const = 0;
461 : /// @}
462 :
463 : /// @{
464 : /*!
465 : * Computes \f$\chi=\partial p / \partial \rho |_{\epsilon}\f$ from the
466 : * \f$\rho\f$ and \f$\epsilon\f$, where \f$p\f$ is the pressure, \f$\rho\f$ is
467 : * the rest mass density, and \f$\epsilon\f$ is the specific internal energy.
468 : */
469 1 : virtual Scalar<double> chi_from_density_and_energy(
470 : const Scalar<double>& /*rest_mass_density*/,
471 : const Scalar<double>& /*specific_internal_energy*/) const = 0;
472 1 : virtual Scalar<DataVector> chi_from_density_and_energy(
473 : const Scalar<DataVector>& /*rest_mass_density*/,
474 : const Scalar<DataVector>& /*specific_internal_energy*/) const = 0;
475 : /// @}
476 :
477 : /// @{
478 : /*!
479 : * Computes \f$\kappa p/\rho^2=(p/\rho^2)\partial p / \partial \epsilon
480 : * |_{\rho}\f$ from \f$\rho\f$ and \f$\epsilon\f$, where \f$p\f$ is the
481 : * pressure, \f$\rho\f$ is the rest mass density, and \f$\epsilon\f$ is the
482 : * specific internal energy.
483 : *
484 : * The reason for not returning just
485 : * \f$\kappa=\partial p / \partial \epsilon\f$ is to avoid division by zero
486 : * for small values of \f$\rho\f$ when assembling the speed of sound with
487 : * some equations of state.
488 : */
489 1 : virtual Scalar<double> kappa_times_p_over_rho_squared_from_density_and_energy(
490 : const Scalar<double>& /*rest_mass_density*/,
491 : const Scalar<double>& /*specific_internal_energy*/) const = 0;
492 : virtual Scalar<DataVector>
493 1 : kappa_times_p_over_rho_squared_from_density_and_energy(
494 : const Scalar<DataVector>& /*rest_mass_density*/,
495 : const Scalar<DataVector>& /*specific_internal_energy*/) const = 0;
496 : /// @}
497 :
498 : /// The lower bound of the electron fraction that is valid for this EOS
499 1 : virtual double electron_fraction_lower_bound() const { return 0.0; }
500 :
501 : /// The upper bound of the electron fraction that is valid for this EOS
502 1 : virtual double electron_fraction_upper_bound() const { return 1.0; }
503 :
504 : /// The lower bound of the rest mass density that is valid for this EOS
505 1 : virtual double rest_mass_density_lower_bound() const = 0;
506 :
507 : /// The upper bound of the rest mass density that is valid for this EOS
508 1 : virtual double rest_mass_density_upper_bound() const = 0;
509 :
510 : /// The lower bound of the specific internal energy that is valid for this EOS
511 : /// at the given rest mass density \f$\rho\f$
512 1 : virtual double specific_internal_energy_lower_bound(
513 : const double rest_mass_density) const = 0;
514 :
515 : /// The upper bound of the specific internal energy that is valid for this EOS
516 : /// at the given rest mass density \f$\rho\f$
517 1 : virtual double specific_internal_energy_upper_bound(
518 : const double rest_mass_density) const = 0;
519 :
520 : /// The lower bound of the temperature that is valid for this EOS
521 1 : virtual double temperature_lower_bound() const { return 0.0; };
522 :
523 : /// The upper bound of the temperature that is valid for this EOS
524 1 : virtual double temperature_upper_bound() const {
525 : return std::numeric_limits<double>::max();
526 : };
527 :
528 : /// The lower bound of the specific enthalpy that is valid for this EOS
529 1 : virtual double specific_enthalpy_lower_bound() const = 0;
530 :
531 : /// The vacuum mass of a baryon for this EOS
532 1 : virtual double baryon_mass() const {
533 : return hydro::units::geometric::default_baryon_mass;
534 : }
535 : };
536 :
537 : /*!
538 : * \ingroup EquationsOfStateGroup
539 : * \brief Base class for equations of state which need three independent
540 : * thermodynamic variables in order to determine the pressure.
541 : *
542 : * The template parameter `IsRelativistic` is `true` for relativistic equations
543 : * of state and `false` for non-relativistic equations of state.
544 : */
545 : template <bool IsRelativistic>
546 1 : class EquationOfState<IsRelativistic, 3> : public PUP::able {
547 : public:
548 0 : static constexpr bool is_relativistic = IsRelativistic;
549 0 : static constexpr size_t thermodynamic_dim = 3;
550 0 : using creatable_classes =
551 : typename detail::DerivedClasses<IsRelativistic, 3>::type;
552 :
553 0 : EquationOfState() = default;
554 0 : EquationOfState(const EquationOfState&) = default;
555 0 : EquationOfState& operator=(const EquationOfState&) = default;
556 0 : EquationOfState(EquationOfState&&) = default;
557 0 : EquationOfState& operator=(EquationOfState&&) = default;
558 0 : ~EquationOfState() override = default;
559 :
560 0 : explicit EquationOfState(CkMigrateMessage* msg) : PUP::able(msg) {}
561 :
562 0 : WRAPPED_PUPable_abstract(EquationOfState); // NOLINT
563 :
564 0 : virtual inline std::unique_ptr<EquationOfState<IsRelativistic, 3>> get_clone()
565 : const = 0;
566 :
567 0 : virtual bool is_equal(
568 : const EquationOfState<IsRelativistic, 3>& rhs) const = 0;
569 :
570 : virtual std::unique_ptr<EquationOfState<IsRelativistic, 3>>
571 0 : promote_to_3d_eos() {
572 : return this->get_clone();
573 : }
574 :
575 : /// \brief Returns `true` if the EOS is barotropic
576 1 : virtual bool is_barotropic() const = 0;
577 :
578 : /// \brief Returns `true` if the EOS is in beta-equilibrium
579 1 : virtual bool is_equilibrium() const = 0;
580 :
581 : /// @{
582 : /*!
583 : * Computes the electron fraction in beta-equilibrium \f$Y_e^{\rm eq}\f$ from
584 : * the rest mass density \f$\rho\f$ and the temperature \f$T\f$.
585 : */
586 1 : virtual Scalar<double> equilibrium_electron_fraction_from_density_temperature(
587 : const Scalar<double>& /*rest_mass_density*/,
588 : const Scalar<double>& /*temperature*/) const = 0;
589 :
590 : virtual Scalar<DataVector>
591 1 : equilibrium_electron_fraction_from_density_temperature(
592 : const Scalar<DataVector>& /*rest_mass_density*/,
593 : const Scalar<DataVector>& /*temperature*/) const = 0;
594 : /// @}
595 :
596 : /// @{
597 : /*!
598 : * Computes the pressure \f$p\f$ from the rest mass density \f$\rho\f$, the
599 : * specific internal energy \f$\epsilon\f$ and electron fraction \f$Y_e\f$.
600 : */
601 1 : virtual Scalar<double> pressure_from_density_and_energy(
602 : const Scalar<double>& /*rest_mass_density*/,
603 : const Scalar<double>& /*specific_internal_energy*/,
604 : const Scalar<double>& /*electron_fraction*/) const = 0;
605 1 : virtual Scalar<DataVector> pressure_from_density_and_energy(
606 : const Scalar<DataVector>& /*rest_mass_density*/,
607 : const Scalar<DataVector>& /*specific_internal_energy*/,
608 : const Scalar<DataVector>& /*electron_fraction*/) const = 0;
609 : /// @}
610 :
611 : /// @{
612 : /*!
613 : * Computes the pressure \f$p\f$ from the rest mass density \f$\rho\f$, the
614 : * temperature \f$T\f$, and electron fraction \f$Y_e\f$.
615 : */
616 1 : virtual Scalar<double> pressure_from_density_and_temperature(
617 : const Scalar<double>& /*rest_mass_density*/,
618 : const Scalar<double>& /*temperature*/,
619 : const Scalar<double>& /*electron_fraction*/) const = 0;
620 1 : virtual Scalar<DataVector> pressure_from_density_and_temperature(
621 : const Scalar<DataVector>& /*rest_mass_density*/,
622 : const Scalar<DataVector>& /*temperature*/,
623 : const Scalar<DataVector>& /*electron_fraction*/) const = 0;
624 : /// @}
625 :
626 : /// @{
627 : /*!
628 : * Computes the temperature \f$T\f$ from the rest mass
629 : * density \f$\rho\f$, the specific internal energy \f$\epsilon\f$,
630 : * and electron fraction \f$Y_e\f$.
631 : */
632 1 : virtual Scalar<double> temperature_from_density_and_energy(
633 : const Scalar<double>& /*rest_mass_density*/,
634 : const Scalar<double>& /*specific_internal_energy*/,
635 : const Scalar<double>& /*electron_fraction*/) const = 0;
636 1 : virtual Scalar<DataVector> temperature_from_density_and_energy(
637 : const Scalar<DataVector>& /*rest_mass_density*/,
638 : const Scalar<DataVector>& /*specific_internal_energy*/,
639 : const Scalar<DataVector>& /*electron_fraction*/) const = 0;
640 : /// @}
641 :
642 : /// @{
643 : /*!
644 : * Computes the specific internal energy \f$\epsilon\f$ from the rest mass
645 : * density \f$\rho\f$, the temperature \f$T\f$, and electron fraction
646 : * \f$Y_e\f$.
647 : */
648 1 : virtual Scalar<double> specific_internal_energy_from_density_and_temperature(
649 : const Scalar<double>& /*rest_mass_density*/,
650 : const Scalar<double>& /*temperature*/,
651 : const Scalar<double>& /*electron_fraction*/
652 : ) const = 0;
653 : virtual Scalar<DataVector>
654 1 : specific_internal_energy_from_density_and_temperature(
655 : const Scalar<DataVector>& /*rest_mass_density*/,
656 : const Scalar<DataVector>& /*temperature*/,
657 : const Scalar<DataVector>& /*electron_fraction*/
658 : ) const = 0;
659 : /// @}
660 :
661 : /// @{
662 : /*!
663 : * Computes adiabatic sound speed squared
664 : * \f[
665 : * c_s^2 \equiv \frac{\partial p}{\partial e} |_{s, Y_e} =
666 : * \frac{\rho}{h}\frac{\partial p}{\partial \rho} |_{e, Y_e} +
667 : * \frac{\partial p}{\partial e}|_{\rho, Y_e}
668 : * \f].
669 : * With \f$p, e\f$ the pressure and energy density respectively,
670 : * \f$s\f$ the entropy density, \f$Y_e\f$ the electron fraction
671 : * \f$\rho\f$ the rest-mass density, and \f$h\f$ the enthalpy density.
672 : * Note that \f$e\f$ is the total energy density and not the internal energy,
673 : * therefore
674 : * \f[
675 : * \frac{\partial p}{\partial \rho} |_{e, Y_e} \neq \chi \equiv \frac{\partial
676 : * p}{\partial \rho} |_{\epsilon, Y_e}
677 : * \f]
678 : * as defined in the 2-d EoS above. By definition
679 : * \f$ e = (1+\epsilon) \rho \f$ so holding \f$e\f$ constant
680 : * \f[
681 : * 0 = \frac{d e}{d \rho} = \frac{\partial e}{\partial \rho} +
682 : * \frac{\partial e}{\partial \epsilon} \frac{\partial \epsilon}{\rho}.
683 : * \f]
684 : * (where we have suppressed \f$ Y_e\f$ dependence)
685 : * So \f$ \partial \epsilon / \partial \rho |_{e} = (1 + \epsilon)/\rho \f$
686 : * and we can expand
687 : * \f[
688 : * \frac{\partial p}{\partial \rho} |_{e, Y_e} = \frac{\partial e}{\partial
689 : * \rho}_{\epsilon, Y_e} + \frac{(1 + \epsilon)}{\rho} \frac{\partial
690 : * e}{\partial \epsilon}|_{\rho, Y_e}
691 : * \f]
692 : * Finally, we can rewrite the entire sound speed using only the rest-mass
693 : * density, specific internal energy, and electron fraction as variables,
694 : * by using \f$ \frac{\partial e}{\partial \epsilon}|_{\rho, Y_e} = 1 \f$
695 : * \f[
696 : * c_s^2 =
697 : * \frac{\rho}{h}\frac{\partial p}{\partial \rho} |_{\epsilon, Y_e} +
698 : * \frac{1}{\rho} \frac{\partial p}{\partial \epsilon}|_{\rho, Y_e} \left(
699 : * 1 - \frac{(1 + \epsilon)\rho}{h}\right)
700 : * \f]
701 : * Which reduces to our preferred form
702 : * \f[
703 : * c_s^2 =
704 : * \frac{\rho}{h}(\chi + \kappa)
705 : * \f]
706 : *
707 : * Computed as a function of temperature, rest-mass density and electron
708 : * fraction. Note that this will break thermodynamic consistency if the
709 : * pressure and internal energy interpolated separately. The precise impact of
710 : * this will depend on the EoS and numerical scheme used for the evolution.
711 : */
712 1 : virtual Scalar<double> sound_speed_squared_from_density_and_temperature(
713 : const Scalar<double>& /*rest_mass_density*/,
714 : const Scalar<double>& /*temperature*/,
715 : const Scalar<double>& /*electron_fraction*/) const = 0;
716 1 : virtual Scalar<DataVector> sound_speed_squared_from_density_and_temperature(
717 : const Scalar<DataVector>& /*rest_mass_density*/,
718 : const Scalar<DataVector>& /*temperature*/,
719 : const Scalar<DataVector>& /*electron_fraction*/) const = 0;
720 : /// @}
721 :
722 : /// The lower bound of the electron fraction that is valid for this EOS
723 1 : virtual double electron_fraction_lower_bound() const = 0;
724 :
725 : /// The upper bound of the electron fraction that is valid for this EOS
726 1 : virtual double electron_fraction_upper_bound() const = 0;
727 :
728 : /// The lower bound of the rest mass density that is valid for this EOS
729 1 : virtual double rest_mass_density_lower_bound() const = 0;
730 :
731 : /// The upper bound of the rest mass density that is valid for this EOS
732 1 : virtual double rest_mass_density_upper_bound() const = 0;
733 :
734 : /// The lower bound of the specific internal energy that is valid for this EOS
735 : /// at the given rest mass density \f$\rho\f$ and electron fraction \f$Y_e\f$.
736 1 : virtual double specific_internal_energy_lower_bound(
737 : const double rest_mass_density, const double electron_fraction) const = 0;
738 :
739 : /// The upper bound of the specific internal energy that is valid for this EOS
740 : /// at the given rest mass density \f$\rho\f$ and electron fraction \f$Y_e\f$.
741 1 : virtual double specific_internal_energy_upper_bound(
742 : const double rest_mass_density, const double electron_fraction) const = 0;
743 :
744 : /// The lower bound of the specific enthalpy that is valid for this EOS
745 1 : virtual double specific_enthalpy_lower_bound() const = 0;
746 :
747 : /// The lower bound of the temperature that is valid for this EOS
748 1 : virtual double temperature_lower_bound() const = 0;
749 :
750 : /// The upper bound of the temperature that is valid for this EOS
751 1 : virtual double temperature_upper_bound() const = 0;
752 :
753 : /// The vacuum mass of a baryon for this EOS
754 1 : virtual double baryon_mass() const {
755 : return hydro::units::geometric::default_baryon_mass;
756 : }
757 : };
758 :
759 : /// Compare two equations of state for equality
760 : template <bool IsRelLhs, bool IsRelRhs, size_t ThermoDimLhs,
761 : size_t ThermoDimRhs>
762 1 : bool operator==(const EquationOfState<IsRelLhs, ThermoDimLhs>& lhs,
763 : const EquationOfState<IsRelRhs, ThermoDimRhs>& rhs) {
764 : if constexpr (IsRelLhs == IsRelRhs and ThermoDimLhs == ThermoDimRhs) {
765 : return typeid(lhs) == typeid(rhs) and lhs.is_equal(rhs);
766 : } else {
767 : return false;
768 : }
769 : }
770 : template <bool IsRelLhs, bool IsRelRhs, size_t ThermoDimLhs,
771 : size_t ThermoDimRhs>
772 0 : bool operator!=(const EquationOfState<IsRelLhs, ThermoDimLhs>& lhs,
773 : const EquationOfState<IsRelRhs, ThermoDimRhs>& rhs) {
774 : return not(lhs == rhs);
775 : }
776 : } // namespace EquationsOfState
777 :
778 : /// \cond
779 : #define EQUATION_OF_STATE_FUNCTIONS_1D \
780 : (pressure_from_density, rest_mass_density_from_enthalpy, \
781 : specific_internal_energy_from_density, chi_from_density, \
782 : kappa_times_p_over_rho_squared_from_density)
783 :
784 : #define EQUATION_OF_STATE_FUNCTIONS_2D \
785 : (pressure_from_density_and_energy, pressure_from_density_and_enthalpy, \
786 : specific_internal_energy_from_density_and_pressure, \
787 : temperature_from_density_and_energy, \
788 : specific_internal_energy_from_density_and_temperature, \
789 : chi_from_density_and_energy, \
790 : kappa_times_p_over_rho_squared_from_density_and_energy)
791 :
792 : #define EQUATION_OF_STATE_FUNCTIONS_3D \
793 : (pressure_from_density_and_energy, pressure_from_density_and_temperature, \
794 : temperature_from_density_and_energy, \
795 : specific_internal_energy_from_density_and_temperature, \
796 : sound_speed_squared_from_density_and_temperature)
797 :
798 : #define EQUATION_OF_STATE_ARGUMENTS_EXPAND(z, n, type) \
799 : BOOST_PP_COMMA_IF(n) const Scalar<type>&
800 :
801 : #define EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS_HELPER(r, DIM, \
802 : FUNCTION_NAME) \
803 : Scalar<double> FUNCTION_NAME(BOOST_PP_REPEAT( \
804 : DIM, EQUATION_OF_STATE_ARGUMENTS_EXPAND, double)) const override; \
805 : Scalar<DataVector> FUNCTION_NAME(BOOST_PP_REPEAT( \
806 : DIM, EQUATION_OF_STATE_ARGUMENTS_EXPAND, DataVector)) const override;
807 :
808 : /// \endcond
809 :
810 : /*!
811 : * \ingroup EquationsOfStateGroup
812 : * \brief Macro used to generate forward declarations of member functions in
813 : * derived classes
814 : */
815 1 : #define EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(DERIVED, DIM) \
816 : BOOST_PP_LIST_FOR_EACH( \
817 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS_HELPER, DIM, \
818 : BOOST_PP_TUPLE_TO_LIST(BOOST_PP_TUPLE_ELEM( \
819 : BOOST_PP_SUB(DIM, 1), \
820 : (EQUATION_OF_STATE_FUNCTIONS_1D, EQUATION_OF_STATE_FUNCTIONS_2D, \
821 : EQUATION_OF_STATE_FUNCTIONS_3D)))) \
822 : \
823 : /* clang-tidy: do not use non-const references */ \
824 : void pup(PUP::er& p) override; /* NOLINT */ \
825 : \
826 : explicit DERIVED(CkMigrateMessage* msg);
827 :
828 : /// \cond
829 : #define EQUATION_OF_STATE_FORWARD_ARGUMENTS(z, n, unused) \
830 : BOOST_PP_COMMA_IF(n) arg##n
831 :
832 : #define EQUATION_OF_STATE_ARGUMENTS_EXPAND_NAMED(z, n, type) \
833 : BOOST_PP_COMMA_IF(n) const Scalar<type>& arg##n
834 :
835 : #define EQUATION_OF_STATE_MEMBER_DEFINITIONS_HELPER( \
836 : TEMPLATE, DERIVED, DATA_TYPE, DIM, FUNCTION_NAME) \
837 : TEMPLATE \
838 : Scalar<DATA_TYPE> DERIVED::FUNCTION_NAME(BOOST_PP_REPEAT( \
839 : DIM, EQUATION_OF_STATE_ARGUMENTS_EXPAND_NAMED, DATA_TYPE)) const { \
840 : return FUNCTION_NAME##_impl( \
841 : BOOST_PP_REPEAT(DIM, EQUATION_OF_STATE_FORWARD_ARGUMENTS, UNUSED)); \
842 : }
843 :
844 : #define EQUATION_OF_STATE_MEMBER_DEFINITIONS_HELPER_2(r, ARGS, FUNCTION_NAME) \
845 : EQUATION_OF_STATE_MEMBER_DEFINITIONS_HELPER( \
846 : BOOST_PP_TUPLE_ELEM(0, ARGS), BOOST_PP_TUPLE_ELEM(1, ARGS), \
847 : BOOST_PP_TUPLE_ELEM(2, ARGS), BOOST_PP_TUPLE_ELEM(3, ARGS), \
848 : FUNCTION_NAME)
849 : /// \endcond
850 :
851 : #define EQUATION_OF_STATE_MEMBER_DEFINITIONS(TEMPLATE, DERIVED, DATA_TYPE, \
852 0 : DIM) \
853 : BOOST_PP_LIST_FOR_EACH( \
854 : EQUATION_OF_STATE_MEMBER_DEFINITIONS_HELPER_2, \
855 : (TEMPLATE, DERIVED, DATA_TYPE, DIM), \
856 : BOOST_PP_TUPLE_TO_LIST(BOOST_PP_TUPLE_ELEM( \
857 : BOOST_PP_SUB(DIM, 1), \
858 : (EQUATION_OF_STATE_FUNCTIONS_1D, EQUATION_OF_STATE_FUNCTIONS_2D, \
859 : EQUATION_OF_STATE_FUNCTIONS_3D))))
860 :
861 : /// \cond
862 : #define EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS_HELPER(r, DIM, \
863 : FUNCTION_NAME) \
864 : template <class DataType> \
865 : Scalar<DataType> FUNCTION_NAME##_impl(BOOST_PP_REPEAT( \
866 : DIM, EQUATION_OF_STATE_ARGUMENTS_EXPAND, DataType)) const;
867 : /// \endcond
868 :
869 0 : #define EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(DIM) \
870 : BOOST_PP_LIST_FOR_EACH( \
871 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS_HELPER, DIM, \
872 : BOOST_PP_TUPLE_TO_LIST(BOOST_PP_TUPLE_ELEM( \
873 : BOOST_PP_SUB(DIM, 1), \
874 : (EQUATION_OF_STATE_FUNCTIONS_1D, EQUATION_OF_STATE_FUNCTIONS_2D, \
875 : EQUATION_OF_STATE_FUNCTIONS_3D))))
|
