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 entropy \f$s\f$ from the rest mass density
227 : * \f$\rho\f$.
228 : */
229 1 : virtual Scalar<double> specific_entropy_from_density(
230 : const Scalar<double>& /*rest_mass_density*/) const {
231 : return Scalar<double>{0.0};
232 : }
233 1 : virtual Scalar<DataVector> specific_entropy_from_density(
234 : const Scalar<DataVector>& rest_mass_density) const {
235 : return make_with_value<Scalar<DataVector>>(rest_mass_density, 0.0);
236 : }
237 : /// @}
238 :
239 : /// @{
240 : /*!
241 : * Computes the specific entropy \f$s\f$ from the specific internal energy
242 : * \f$\epsilon\f$.
243 : */
244 1 : virtual Scalar<double> specific_entropy_from_specific_internal_energy(
245 : const Scalar<double>& /*specific_internal_energy*/) const {
246 : return Scalar<double>{0.0};
247 : }
248 1 : virtual Scalar<DataVector> specific_entropy_from_specific_internal_energy(
249 : const Scalar<DataVector>& specific_internal_energy) const {
250 : return make_with_value<Scalar<DataVector>>(specific_internal_energy, 0.0);
251 : }
252 : /// @}
253 :
254 : /// @{
255 : /*!
256 : * Computes the specific internal energy \f$\epsilon\f$ from the rest mass
257 : * density \f$\rho\f$.
258 : */
259 1 : virtual Scalar<double> specific_internal_energy_from_density(
260 : const Scalar<double>& /*rest_mass_density*/) const = 0;
261 1 : virtual Scalar<DataVector> specific_internal_energy_from_density(
262 : const Scalar<DataVector>& /*rest_mass_density*/) const = 0;
263 : /// @}
264 :
265 : /// @{
266 : /*!
267 : * Computes the temperature \f$T\f$ from the rest mass
268 : * density \f$\rho\f$.
269 : */
270 1 : virtual Scalar<double> temperature_from_density(
271 : const Scalar<double>& /*rest_mass_density*/) const {
272 : return Scalar<double>{0.0};
273 : }
274 1 : virtual Scalar<DataVector> temperature_from_density(
275 : const Scalar<DataVector>& rest_mass_density) const {
276 : return make_with_value<Scalar<DataVector>>(rest_mass_density, 0.0);
277 : }
278 : /// @}
279 :
280 : /// @{
281 : /*!
282 : * Computes the temperature \f$\T\f$ from the specific internal energy
283 : * \f$\epsilon\f$.
284 : */
285 1 : virtual Scalar<double> temperature_from_specific_internal_energy(
286 : const Scalar<double>& /*specific_internal_energy*/) const {
287 : return Scalar<double>{0.0};
288 : }
289 1 : virtual Scalar<DataVector> temperature_from_specific_internal_energy(
290 : const Scalar<DataVector>& specific_internal_energy) const {
291 : return make_with_value<Scalar<DataVector>>(specific_internal_energy, 0.0);
292 : }
293 : /// @}
294 :
295 : /// @{
296 : /*!
297 : * Computes \f$\chi=\partial p / \partial \rho\f$ from \f$\rho\f$, where
298 : * \f$p\f$ is the pressure and \f$\rho\f$ is the rest mass density.
299 : */
300 1 : virtual Scalar<double> chi_from_density(
301 : const Scalar<double>& /*rest_mass_density*/) const = 0;
302 1 : virtual Scalar<DataVector> chi_from_density(
303 : const Scalar<DataVector>& /*rest_mass_density*/) const = 0;
304 : /// @}
305 :
306 : /// @{
307 : /*!
308 : * Computes \f$\kappa p/\rho^2=(p/\rho^2)\partial p / \partial \epsilon\f$
309 : * from \f$\rho\f$, where \f$p\f$ is the pressure, \f$\rho\f$ is the rest mass
310 : * density, and \f$\epsilon\f$ is the specific internal energy.
311 : *
312 : * The reason for not returning just
313 : * \f$\kappa=\partial p / \partial \epsilon\f$ is to avoid division by zero
314 : * for small values of \f$\rho\f$ when assembling the speed of sound with
315 : * some equations of state.
316 : */
317 1 : virtual Scalar<double> kappa_times_p_over_rho_squared_from_density(
318 : const Scalar<double>& /*rest_mass_density*/) const = 0;
319 1 : virtual Scalar<DataVector> kappa_times_p_over_rho_squared_from_density(
320 : const Scalar<DataVector>& /*rest_mass_density*/) const = 0;
321 :
322 : /// The lower bound of the electron fraction that is valid for this EOS
323 1 : virtual double electron_fraction_lower_bound() const { return 0.0; }
324 :
325 : /// The upper bound of the electron fraction that is valid for this EOS
326 1 : virtual double electron_fraction_upper_bound() const { return 1.0; }
327 :
328 : /// The lower bound of the rest mass density that is valid for this EOS
329 1 : virtual double rest_mass_density_lower_bound() const = 0;
330 :
331 : /// The upper bound of the rest mass density that is valid for this EOS
332 1 : virtual double rest_mass_density_upper_bound() const = 0;
333 :
334 : /// The lower bound of the temperature that is valid for this EOS
335 1 : virtual double temperature_lower_bound() const { return 0.0; };
336 :
337 : /// The upper bound of the temperature that is valid for this EOS
338 1 : virtual double temperature_upper_bound() const {
339 : return std::numeric_limits<double>::max();
340 : };
341 :
342 : /// The lower bound of the specific internal energy that is valid for this EOS
343 1 : virtual double specific_internal_energy_lower_bound() const { return 0.0; };
344 :
345 : /// The upper bound of the specific internal energy that is valid for this EOS
346 1 : virtual double specific_internal_energy_upper_bound() const {
347 : return std::numeric_limits<double>::max();
348 : };
349 :
350 : /// The lower bound of the specific enthalpy that is valid for this EOS
351 1 : virtual double specific_enthalpy_lower_bound() const = 0;
352 :
353 : /// The vacuum mass of a baryon for this EOS
354 1 : virtual double baryon_mass() const {
355 : return hydro::units::geometric::default_baryon_mass;
356 : }
357 : };
358 :
359 : /*!
360 : * \ingroup EquationsOfStateGroup
361 : * \brief Base class for equations of state which need two independent
362 : * thermodynamic variables in order to determine the pressure.
363 : *
364 : * The template parameter `IsRelativistic` is `true` for relativistic equations
365 : * of state and `false` for non-relativistic equations of state.
366 : */
367 : template <bool IsRelativistic>
368 1 : class EquationOfState<IsRelativistic, 2> : public PUP::able {
369 : public:
370 0 : static constexpr bool is_relativistic = IsRelativistic;
371 0 : static constexpr size_t thermodynamic_dim = 2;
372 0 : using creatable_classes =
373 : typename detail::DerivedClasses<IsRelativistic, 2>::type;
374 :
375 0 : EquationOfState() = default;
376 0 : EquationOfState(const EquationOfState&) = default;
377 0 : EquationOfState& operator=(const EquationOfState&) = default;
378 0 : EquationOfState(EquationOfState&&) = default;
379 0 : EquationOfState& operator=(EquationOfState&&) = default;
380 0 : ~EquationOfState() override = default;
381 :
382 0 : explicit EquationOfState(CkMigrateMessage* msg) : PUP::able(msg) {}
383 :
384 0 : WRAPPED_PUPable_abstract(EquationOfState); // NOLINT
385 :
386 0 : virtual inline std::unique_ptr<EquationOfState<IsRelativistic, 2>> get_clone()
387 : const = 0;
388 :
389 0 : virtual bool is_equal(
390 : const EquationOfState<IsRelativistic, 2>& rhs) const = 0;
391 :
392 : virtual std::unique_ptr<EquationOfState<IsRelativistic, 2>>
393 0 : promote_to_2d_eos() const {
394 : return this->get_clone();
395 : }
396 :
397 : virtual std::unique_ptr<EquationOfState<IsRelativistic, 3>>
398 0 : promote_to_3d_eos() const = 0;
399 :
400 : /// \brief Returns `true` if the EOS is barotropic
401 1 : virtual bool is_barotropic() const = 0;
402 :
403 : /// \brief Returns `true` if the EOS is in beta-equilibrium
404 1 : virtual bool is_equilibrium() const {
405 : return true;
406 : }
407 :
408 : /// @{
409 : /*!
410 : * Computes the electron fraction in beta-equilibrium \f$Y_e^{\rm eq}\f$ from
411 : * the rest mass density \f$\rho\f$ and the temperature \f$T\f$.
412 : */
413 1 : virtual Scalar<double> equilibrium_electron_fraction_from_density_temperature(
414 : const Scalar<double>& rest_mass_density,
415 : const Scalar<double>& /*temperature*/) const {
416 : return make_with_value<Scalar<double>>(rest_mass_density, 0.1);
417 : }
418 :
419 : virtual Scalar<DataVector>
420 1 : equilibrium_electron_fraction_from_density_temperature(
421 : const Scalar<DataVector>& rest_mass_density,
422 : const Scalar<DataVector>& /*temperature*/) const {
423 : return make_with_value<Scalar<DataVector>>(rest_mass_density, 0.1);
424 : }
425 : /// @}
426 :
427 : /// @{
428 : /*!
429 : * Computes the pressure \f$p\f$ from the rest mass density \f$\rho\f$ and the
430 : * specific internal energy \f$\epsilon\f$.
431 : */
432 1 : virtual Scalar<double> pressure_from_density_and_energy(
433 : const Scalar<double>& /*rest_mass_density*/,
434 : const Scalar<double>& /*specific_internal_energy*/) const = 0;
435 1 : virtual Scalar<DataVector> pressure_from_density_and_energy(
436 : const Scalar<DataVector>& /*rest_mass_density*/,
437 : const Scalar<DataVector>& /*specific_internal_energy*/) const = 0;
438 : /// @}
439 :
440 : /// @{
441 : /*!
442 : * Computes the pressure \f$p\f$ from the rest mass density \f$\rho\f$ and the
443 : * specific enthalpy \f$h\f$.
444 : */
445 1 : virtual Scalar<double> pressure_from_density_and_enthalpy(
446 : const Scalar<double>& /*rest_mass_density*/,
447 : const Scalar<double>& /*specific_enthalpy*/) const = 0;
448 1 : virtual Scalar<DataVector> pressure_from_density_and_enthalpy(
449 : const Scalar<DataVector>& /*rest_mass_density*/,
450 : const Scalar<DataVector>& /*specific_enthalpy*/) const = 0;
451 : /// @}
452 :
453 : /// @{
454 : /*!
455 : * Computes the specific entropy \f$s\f$ from the rest mass density \f$\rho\f$
456 : * and the specific internal energy \f$\epsilon\f$.
457 : */
458 1 : virtual Scalar<double> specific_entropy_from_density_and_energy(
459 : const Scalar<double>& /*rest_mass_density*/,
460 : const Scalar<double>& /*specific_internal_energy*/) const = 0;
461 :
462 1 : virtual Scalar<DataVector> specific_entropy_from_density_and_energy(
463 : const Scalar<DataVector>& /*rest_mass_density*/,
464 : const Scalar<DataVector>& /*specific_internal_energy*/) const = 0;
465 : /// @}
466 :
467 : /// @{
468 : /*!
469 : * Computes the specific entropy \f$s\f$ from the rest mass density \f$\rho\f$
470 : * and the temperature \f$T\f$.
471 : */
472 1 : virtual Scalar<double> specific_entropy_from_density_and_temperature(
473 : const Scalar<double>& /*rest_mass_density*/,
474 : const Scalar<double>& /*temperature*/) const = 0;
475 :
476 1 : virtual Scalar<DataVector> specific_entropy_from_density_and_temperature(
477 : const Scalar<DataVector>& /*rest_mass_density*/,
478 : const Scalar<DataVector>& /*temperature*/) const = 0;
479 : /// @}
480 :
481 : /// @{
482 : /*!
483 : * Computes the specific internal energy \f$\epsilon\f$ from the rest mass
484 : * density \f$\rho\f$ and the pressure \f$p\f$.
485 : */
486 1 : virtual Scalar<double> specific_internal_energy_from_density_and_pressure(
487 : const Scalar<double>& /*rest_mass_density*/,
488 : const Scalar<double>& /*pressure*/) const = 0;
489 1 : virtual Scalar<DataVector> specific_internal_energy_from_density_and_pressure(
490 : const Scalar<DataVector>& /*rest_mass_density*/,
491 : const Scalar<DataVector>& /*pressure*/) const = 0;
492 : /// @}
493 :
494 : /// @{
495 : /*!
496 : * Computes the temperature \f$T\f$ from the rest mass
497 : * density \f$\rho\f$ and the specific internal energy \f$\epsilon\f$.
498 : */
499 1 : virtual Scalar<double> temperature_from_density_and_energy(
500 : const Scalar<double>& /*rest_mass_density*/,
501 : const Scalar<double>& /*specific_internal_energy*/) const = 0;
502 1 : virtual Scalar<DataVector> temperature_from_density_and_energy(
503 : const Scalar<DataVector>& /*rest_mass_density*/,
504 : const Scalar<DataVector>& /*specific_internal_energy*/) const = 0;
505 : /// @}
506 :
507 : /// @{
508 : /*!
509 : * Computes the specific internal energy \f$\epsilon\f$ from the rest mass
510 : * density \f$\rho\f$ and the temperature \f$T\f$.
511 : */
512 1 : virtual Scalar<double> specific_internal_energy_from_density_and_temperature(
513 : const Scalar<double>& /*rest_mass_density*/,
514 : const Scalar<double>& /*temperature*/) const = 0;
515 : virtual Scalar<DataVector>
516 1 : specific_internal_energy_from_density_and_temperature(
517 : const Scalar<DataVector>& /*rest_mass_density*/,
518 : const Scalar<DataVector>& /*temperature*/) const = 0;
519 : /// @}
520 :
521 : /// @{
522 : /*!
523 : * Computes \f$\chi=\partial p / \partial \rho |_{\epsilon}\f$ from the
524 : * \f$\rho\f$ and \f$\epsilon\f$, where \f$p\f$ is the pressure, \f$\rho\f$ is
525 : * the rest mass density, and \f$\epsilon\f$ is the specific internal energy.
526 : */
527 1 : virtual Scalar<double> chi_from_density_and_energy(
528 : const Scalar<double>& /*rest_mass_density*/,
529 : const Scalar<double>& /*specific_internal_energy*/) const = 0;
530 1 : virtual Scalar<DataVector> chi_from_density_and_energy(
531 : const Scalar<DataVector>& /*rest_mass_density*/,
532 : const Scalar<DataVector>& /*specific_internal_energy*/) const = 0;
533 : /// @}
534 :
535 : /// @{
536 : /*!
537 : * Computes \f$\kappa p/\rho^2=(p/\rho^2)\partial p / \partial \epsilon
538 : * |_{\rho}\f$ from \f$\rho\f$ and \f$\epsilon\f$, where \f$p\f$ is the
539 : * pressure, \f$\rho\f$ is the rest mass density, and \f$\epsilon\f$ is the
540 : * specific internal energy.
541 : *
542 : * The reason for not returning just
543 : * \f$\kappa=\partial p / \partial \epsilon\f$ is to avoid division by zero
544 : * for small values of \f$\rho\f$ when assembling the speed of sound with
545 : * some equations of state.
546 : */
547 1 : virtual Scalar<double> kappa_times_p_over_rho_squared_from_density_and_energy(
548 : const Scalar<double>& /*rest_mass_density*/,
549 : const Scalar<double>& /*specific_internal_energy*/) const = 0;
550 : virtual Scalar<DataVector>
551 1 : kappa_times_p_over_rho_squared_from_density_and_energy(
552 : const Scalar<DataVector>& /*rest_mass_density*/,
553 : const Scalar<DataVector>& /*specific_internal_energy*/) const = 0;
554 : /// @}
555 :
556 : /// The lower bound of the electron fraction that is valid for this EOS
557 1 : virtual double electron_fraction_lower_bound() const { return 0.0; }
558 :
559 : /// The upper bound of the electron fraction that is valid for this EOS
560 1 : virtual double electron_fraction_upper_bound() const { return 1.0; }
561 :
562 : /// The lower bound of the rest mass density that is valid for this EOS
563 1 : virtual double rest_mass_density_lower_bound() const = 0;
564 :
565 : /// The upper bound of the rest mass density that is valid for this EOS
566 1 : virtual double rest_mass_density_upper_bound() const = 0;
567 :
568 : /// The lower bound of the specific internal energy that is valid for this EOS
569 : /// at the given rest mass density \f$\rho\f$
570 1 : virtual double specific_internal_energy_lower_bound(
571 : const double rest_mass_density) const = 0;
572 :
573 : /// The upper bound of the specific internal energy that is valid for this EOS
574 : /// at the given rest mass density \f$\rho\f$
575 1 : virtual double specific_internal_energy_upper_bound(
576 : const double rest_mass_density) const = 0;
577 :
578 : /// The lower bound of the temperature that is valid for this EOS
579 1 : virtual double temperature_lower_bound() const { return 0.0; };
580 :
581 : /// The upper bound of the temperature that is valid for this EOS
582 1 : virtual double temperature_upper_bound() const {
583 : return std::numeric_limits<double>::max();
584 : };
585 :
586 : /// The lower bound of the specific enthalpy that is valid for this EOS
587 1 : virtual double specific_enthalpy_lower_bound() const = 0;
588 :
589 : /// The vacuum mass of a baryon for this EOS
590 1 : virtual double baryon_mass() const {
591 : return hydro::units::geometric::default_baryon_mass;
592 : }
593 : };
594 :
595 : /*!
596 : * \ingroup EquationsOfStateGroup
597 : * \brief Base class for equations of state which need three independent
598 : * thermodynamic variables in order to determine the pressure.
599 : *
600 : * The template parameter `IsRelativistic` is `true` for relativistic equations
601 : * of state and `false` for non-relativistic equations of state.
602 : */
603 : template <bool IsRelativistic>
604 : class EquationOfState<IsRelativistic, 3> : public PUP::able {
605 : public:
606 : static constexpr bool is_relativistic = IsRelativistic;
607 : static constexpr size_t thermodynamic_dim = 3;
608 : using creatable_classes =
609 : typename detail::DerivedClasses<IsRelativistic, 3>::type;
610 :
611 : EquationOfState() = default;
612 : EquationOfState(const EquationOfState&) = default;
613 : EquationOfState& operator=(const EquationOfState&) = default;
614 : EquationOfState(EquationOfState&&) = default;
615 : EquationOfState& operator=(EquationOfState&&) = default;
616 : ~EquationOfState() override = default;
617 :
618 : explicit EquationOfState(CkMigrateMessage* msg) : PUP::able(msg) {}
619 :
620 : WRAPPED_PUPable_abstract(EquationOfState); // NOLINT
621 :
622 : virtual inline std::unique_ptr<EquationOfState<IsRelativistic, 3>> get_clone()
623 : const = 0;
624 :
625 : virtual bool is_equal(
626 : const EquationOfState<IsRelativistic, 3>& rhs) const = 0;
627 :
628 : virtual std::unique_ptr<EquationOfState<IsRelativistic, 3>>
629 : promote_to_3d_eos() const {
630 : return this->get_clone();
631 : }
632 :
633 : /// \brief Returns `true` if the EOS is barotropic
634 : virtual bool is_barotropic() const = 0;
635 :
636 : /// \brief Returns `true` if the EOS is in beta-equilibrium
637 : virtual bool is_equilibrium() const = 0;
638 :
639 : /// @{
640 : /*!
641 : * Computes the electron fraction in beta-equilibrium \f$Y_e^{\rm eq}\f$ from
642 : * the rest mass density \f$\rho\f$ and the temperature \f$T\f$.
643 : */
644 : virtual Scalar<double> equilibrium_electron_fraction_from_density_temperature(
645 : const Scalar<double>& /*rest_mass_density*/,
646 : const Scalar<double>& /*temperature*/) const = 0;
647 :
648 : virtual Scalar<DataVector>
649 : equilibrium_electron_fraction_from_density_temperature(
650 : const Scalar<DataVector>& /*rest_mass_density*/,
651 : const Scalar<DataVector>& /*temperature*/) const = 0;
652 : /// @}
653 :
654 : /// @{
655 : /*!
656 : * Computes the pressure \f$p\f$ from the rest mass density \f$\rho\f$, the
657 : * specific internal energy \f$\epsilon\f$ and electron fraction \f$Y_e\f$.
658 : */
659 : virtual Scalar<double> pressure_from_density_and_energy(
660 : const Scalar<double>& /*rest_mass_density*/,
661 : const Scalar<double>& /*specific_internal_energy*/,
662 : const Scalar<double>& /*electron_fraction*/) const = 0;
663 : virtual Scalar<DataVector> pressure_from_density_and_energy(
664 : const Scalar<DataVector>& /*rest_mass_density*/,
665 : const Scalar<DataVector>& /*specific_internal_energy*/,
666 : const Scalar<DataVector>& /*electron_fraction*/) const = 0;
667 : /// @}
668 :
669 : /// @{
670 : /*!
671 : * Computes the pressure \f$p\f$ from the rest mass density \f$\rho\f$, the
672 : * temperature \f$T\f$, and electron fraction \f$Y_e\f$.
673 : */
674 : virtual Scalar<double> pressure_from_density_and_temperature(
675 : const Scalar<double>& /*rest_mass_density*/,
676 : const Scalar<double>& /*temperature*/,
677 : const Scalar<double>& /*electron_fraction*/) const = 0;
678 : virtual Scalar<DataVector> pressure_from_density_and_temperature(
679 : const Scalar<DataVector>& /*rest_mass_density*/,
680 : const Scalar<DataVector>& /*temperature*/,
681 : const Scalar<DataVector>& /*electron_fraction*/) const = 0;
682 : /// @}
683 :
684 : /// @{
685 : /*!
686 : * Computes the specific entropy \f$s\f$ from the rest mass density \f$\rho\f$
687 : * and the specific internal energy \f$\epsilon\f$ and electron fraction
688 : * \f$Y_e\f$.
689 : */
690 : virtual Scalar<double> specific_entropy_from_density_and_energy(
691 : const Scalar<double>& /*rest_mass_density*/,
692 : const Scalar<double>& /*specific_internal_energy*/,
693 : const Scalar<double>& /*electron_fraction*/) const = 0;
694 :
695 : virtual Scalar<DataVector> specific_entropy_from_density_and_energy(
696 : const Scalar<DataVector>& /*rest_mass_density*/,
697 : const Scalar<DataVector>& /*specific_internal_energy*/,
698 : const Scalar<DataVector>& /*electron_fraction*/) const = 0;
699 : /// @}
700 :
701 : /// @{
702 : /*!
703 : * Computes the specific entropy \f$s\f$ from the rest mass density \f$\rho\f$
704 : * and the temperature \f$T\f$ and electron fraction \f$Y_e\f$.
705 : */
706 : virtual Scalar<double> specific_entropy_from_density_and_temperature(
707 : const Scalar<double>& /*rest_mass_density*/,
708 : const Scalar<double>& /*temperature*/,
709 : const Scalar<double>& /*electron_fraction*/) const = 0;
710 :
711 : virtual Scalar<DataVector> specific_entropy_from_density_and_temperature(
712 : const Scalar<DataVector>& /*rest_mass_density*/,
713 : const Scalar<DataVector>& /*temperature*/,
714 : const Scalar<DataVector>& /*electron_fraction*/) const = 0;
715 : /// @}
716 :
717 : /// @{
718 : /*!
719 : * Computes the temperature \f$T\f$ from the rest mass
720 : * density \f$\rho\f$, the specific internal energy \f$\epsilon\f$,
721 : * and electron fraction \f$Y_e\f$.
722 : */
723 : virtual Scalar<double> temperature_from_density_and_energy(
724 : const Scalar<double>& /*rest_mass_density*/,
725 : const Scalar<double>& /*specific_internal_energy*/,
726 : const Scalar<double>& /*electron_fraction*/) const = 0;
727 : virtual Scalar<DataVector> temperature_from_density_and_energy(
728 : const Scalar<DataVector>& /*rest_mass_density*/,
729 : const Scalar<DataVector>& /*specific_internal_energy*/,
730 : const Scalar<DataVector>& /*electron_fraction*/) const = 0;
731 : /// @}
732 :
733 : /// @{
734 : /*!
735 : * Computes the specific internal energy \f$\epsilon\f$ from the rest mass
736 : * density \f$\rho\f$, the temperature \f$T\f$, and electron fraction
737 : * \f$Y_e\f$.
738 : */
739 : virtual Scalar<double> specific_internal_energy_from_density_and_temperature(
740 : const Scalar<double>& /*rest_mass_density*/,
741 : const Scalar<double>& /*temperature*/,
742 : const Scalar<double>& /*electron_fraction*/
743 : ) const = 0;
744 : virtual Scalar<DataVector>
745 : specific_internal_energy_from_density_and_temperature(
746 : const Scalar<DataVector>& /*rest_mass_density*/,
747 : const Scalar<DataVector>& /*temperature*/,
748 : const Scalar<DataVector>& /*electron_fraction*/
749 : ) const = 0;
750 : /// @}
751 :
752 : /// @{
753 : /*!
754 : * Computes adiabatic sound speed squared
755 : * \f[
756 : * c_s^2 \equiv \frac{\partial p}{\partial e} |_{s, Y_e} =
757 : * \frac{\rho}{h}\frac{\partial p}{\partial \rho} |_{e, Y_e} +
758 : * \frac{\partial p}{\partial e}|_{\rho, Y_e}
759 : * \f].
760 : * With \f$p, e\f$ the pressure and energy density respectively,
761 : * \f$s\f$ the entropy density, \f$Y_e\f$ the electron fraction
762 : * \f$\rho\f$ the rest-mass density, and \f$h\f$ the enthalpy density.
763 : * Note that \f$e\f$ is the total energy density and not the internal energy,
764 : * therefore
765 : * \f[
766 : * \frac{\partial p}{\partial \rho} |_{e, Y_e} \neq \chi \equiv \frac{\partial
767 : * p}{\partial \rho} |_{\epsilon, Y_e}
768 : * \f]
769 : * as defined in the 2-d EoS above. By definition
770 : * \f$ e = (1+\epsilon) \rho \f$ so holding \f$e\f$ constant
771 : * \f[
772 : * 0 = \frac{d e}{d \rho} = \frac{\partial e}{\partial \rho} +
773 : * \frac{\partial e}{\partial \epsilon} \frac{\partial \epsilon}{\rho}.
774 : * \f]
775 : * (where we have suppressed \f$ Y_e\f$ dependence)
776 : * So \f$ \partial \epsilon / \partial \rho |_{e} = (1 + \epsilon)/\rho \f$
777 : * and we can expand
778 : * \f[
779 : * \frac{\partial p}{\partial \rho} |_{e, Y_e} = \frac{\partial e}{\partial
780 : * \rho}_{\epsilon, Y_e} + \frac{(1 + \epsilon)}{\rho} \frac{\partial
781 : * e}{\partial \epsilon}|_{\rho, Y_e}
782 : * \f]
783 : * Finally, we can rewrite the entire sound speed using only the rest-mass
784 : * density, specific internal energy, and electron fraction as variables,
785 : * by using \f$ \frac{\partial e}{\partial \epsilon}|_{\rho, Y_e} = 1 \f$
786 : * \f[
787 : * c_s^2 =
788 : * \frac{\rho}{h}\frac{\partial p}{\partial \rho} |_{\epsilon, Y_e} +
789 : * \frac{1}{\rho} \frac{\partial p}{\partial \epsilon}|_{\rho, Y_e} \left(
790 : * 1 - \frac{(1 + \epsilon)\rho}{h}\right)
791 : * \f]
792 : * Which reduces to our preferred form
793 : * \f[
794 : * c_s^2 =
795 : * \frac{\rho}{h}(\chi + \kappa)
796 : * \f]
797 : *
798 : * Computed as a function of temperature, rest-mass density and electron
799 : * fraction. Note that this will break thermodynamic consistency if the
800 : * pressure and internal energy interpolated separately. The precise impact of
801 : * this will depend on the EoS and numerical scheme used for the evolution.
802 : */
803 : virtual Scalar<double> sound_speed_squared_from_density_and_temperature(
804 : const Scalar<double>& /*rest_mass_density*/,
805 : const Scalar<double>& /*temperature*/,
806 : const Scalar<double>& /*electron_fraction*/) const = 0;
807 : virtual Scalar<DataVector> sound_speed_squared_from_density_and_temperature(
808 : const Scalar<DataVector>& /*rest_mass_density*/,
809 : const Scalar<DataVector>& /*temperature*/,
810 : const Scalar<DataVector>& /*electron_fraction*/) const = 0;
811 : /// @}
812 :
813 : /// The lower bound of the electron fraction that is valid for this EOS
814 : virtual double electron_fraction_lower_bound() const = 0;
815 :
816 : /// The upper bound of the electron fraction that is valid for this EOS
817 : virtual double electron_fraction_upper_bound() const = 0;
818 :
819 : /// The lower bound of the rest mass density that is valid for this EOS
820 : virtual double rest_mass_density_lower_bound() const = 0;
821 :
822 : /// The upper bound of the rest mass density that is valid for this EOS
823 : virtual double rest_mass_density_upper_bound() const = 0;
824 :
825 : /// The lower bound of the specific internal energy that is valid for this EOS
826 : /// at the given rest mass density \f$\rho\f$ and electron fraction \f$Y_e\f$.
827 : virtual double specific_internal_energy_lower_bound(
828 : const double rest_mass_density, const double electron_fraction) const = 0;
829 :
830 : /// The upper bound of the specific internal energy that is valid for this EOS
831 : /// at the given rest mass density \f$\rho\f$ and electron fraction \f$Y_e\f$.
832 : virtual double specific_internal_energy_upper_bound(
833 : const double rest_mass_density, const double electron_fraction) const = 0;
834 :
835 : /// The lower bound of the specific enthalpy that is valid for this EOS
836 : virtual double specific_enthalpy_lower_bound() const = 0;
837 :
838 : /// The lower bound of the temperature that is valid for this EOS
839 : virtual double temperature_lower_bound() const = 0;
840 :
841 : /// The upper bound of the temperature that is valid for this EOS
842 : virtual double temperature_upper_bound() const = 0;
843 :
844 : /// The vacuum mass of a baryon for this EOS
845 : virtual double baryon_mass() const {
846 : return hydro::units::geometric::default_baryon_mass;
847 : }
848 : };
849 :
850 : /// Compare two equations of state for equality
851 : template <bool IsRelLhs, bool IsRelRhs, size_t ThermoDimLhs,
852 : size_t ThermoDimRhs>
853 1 : bool operator==(const EquationOfState<IsRelLhs, ThermoDimLhs>& lhs,
854 : const EquationOfState<IsRelRhs, ThermoDimRhs>& rhs) {
855 : if constexpr (IsRelLhs == IsRelRhs and ThermoDimLhs == ThermoDimRhs) {
856 : return typeid(lhs) == typeid(rhs) and lhs.is_equal(rhs);
857 : } else {
858 : return false;
859 : }
860 : }
861 : template <bool IsRelLhs, bool IsRelRhs, size_t ThermoDimLhs,
862 : size_t ThermoDimRhs>
863 0 : bool operator!=(const EquationOfState<IsRelLhs, ThermoDimLhs>& lhs,
864 : const EquationOfState<IsRelRhs, ThermoDimRhs>& rhs) {
865 : return not(lhs == rhs);
866 : }
867 : } // namespace EquationsOfState
868 :
869 : /// \cond
870 : #define EQUATION_OF_STATE_FUNCTIONS_1D \
871 : (pressure_from_density, rest_mass_density_from_enthalpy, \
872 : specific_internal_energy_from_density, chi_from_density, \
873 : kappa_times_p_over_rho_squared_from_density)
874 :
875 : #define EQUATION_OF_STATE_FUNCTIONS_2D \
876 : (pressure_from_density_and_energy, pressure_from_density_and_enthalpy, \
877 : specific_entropy_from_density_and_energy, \
878 : specific_entropy_from_density_and_temperature, \
879 : specific_internal_energy_from_density_and_pressure, \
880 : temperature_from_density_and_energy, \
881 : specific_internal_energy_from_density_and_temperature, \
882 : chi_from_density_and_energy, \
883 : kappa_times_p_over_rho_squared_from_density_and_energy)
884 :
885 : #define EQUATION_OF_STATE_FUNCTIONS_3D \
886 : (pressure_from_density_and_energy, pressure_from_density_and_temperature, \
887 : specific_entropy_from_density_and_energy, \
888 : specific_entropy_from_density_and_temperature, \
889 : temperature_from_density_and_energy, \
890 : specific_internal_energy_from_density_and_temperature, \
891 : sound_speed_squared_from_density_and_temperature)
892 :
893 : #define EQUATION_OF_STATE_ARGUMENTS_EXPAND(z, n, type) \
894 : BOOST_PP_COMMA_IF(n) const Scalar<type>&
895 :
896 : #define EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS_HELPER(r, DIM, \
897 : FUNCTION_NAME) \
898 : Scalar<double> FUNCTION_NAME(BOOST_PP_REPEAT( \
899 : DIM, EQUATION_OF_STATE_ARGUMENTS_EXPAND, double)) const override; \
900 : Scalar<DataVector> FUNCTION_NAME(BOOST_PP_REPEAT( \
901 : DIM, EQUATION_OF_STATE_ARGUMENTS_EXPAND, DataVector)) const override;
902 :
903 : /// \endcond
904 :
905 : /*!
906 : * \ingroup EquationsOfStateGroup
907 : * \brief Macro used to generate forward declarations of member functions in
908 : * derived classes
909 : */
910 1 : #define EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(DERIVED, DIM) \
911 : BOOST_PP_LIST_FOR_EACH( \
912 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS_HELPER, DIM, \
913 : BOOST_PP_TUPLE_TO_LIST(BOOST_PP_TUPLE_ELEM( \
914 : BOOST_PP_SUB(DIM, 1), \
915 : (EQUATION_OF_STATE_FUNCTIONS_1D, EQUATION_OF_STATE_FUNCTIONS_2D, \
916 : EQUATION_OF_STATE_FUNCTIONS_3D)))) \
917 : \
918 : /* clang-tidy: do not use non-const references */ \
919 : void pup(PUP::er& p) override; /* NOLINT */ \
920 : \
921 : explicit DERIVED(CkMigrateMessage* msg);
922 :
923 : /// \cond
924 : #define EQUATION_OF_STATE_FORWARD_ARGUMENTS(z, n, unused) \
925 : BOOST_PP_COMMA_IF(n) arg##n
926 :
927 : #define EQUATION_OF_STATE_ARGUMENTS_EXPAND_NAMED(z, n, type) \
928 : BOOST_PP_COMMA_IF(n) const Scalar<type>& arg##n
929 :
930 : #define EQUATION_OF_STATE_MEMBER_DEFINITIONS_HELPER( \
931 : TEMPLATE, DERIVED, DATA_TYPE, DIM, FUNCTION_NAME) \
932 : TEMPLATE \
933 : Scalar<DATA_TYPE> DERIVED::FUNCTION_NAME(BOOST_PP_REPEAT( \
934 : DIM, EQUATION_OF_STATE_ARGUMENTS_EXPAND_NAMED, DATA_TYPE)) const { \
935 : return FUNCTION_NAME##_impl( \
936 : BOOST_PP_REPEAT(DIM, EQUATION_OF_STATE_FORWARD_ARGUMENTS, UNUSED)); \
937 : }
938 :
939 : #define EQUATION_OF_STATE_MEMBER_DEFINITIONS_HELPER_2(r, ARGS, FUNCTION_NAME) \
940 : EQUATION_OF_STATE_MEMBER_DEFINITIONS_HELPER( \
941 : BOOST_PP_TUPLE_ELEM(0, ARGS), BOOST_PP_TUPLE_ELEM(1, ARGS), \
942 : BOOST_PP_TUPLE_ELEM(2, ARGS), BOOST_PP_TUPLE_ELEM(3, ARGS), \
943 : FUNCTION_NAME)
944 : /// \endcond
945 :
946 : #define EQUATION_OF_STATE_MEMBER_DEFINITIONS(TEMPLATE, DERIVED, DATA_TYPE, \
947 0 : DIM) \
948 : BOOST_PP_LIST_FOR_EACH( \
949 : EQUATION_OF_STATE_MEMBER_DEFINITIONS_HELPER_2, \
950 : (TEMPLATE, DERIVED, DATA_TYPE, DIM), \
951 : BOOST_PP_TUPLE_TO_LIST(BOOST_PP_TUPLE_ELEM( \
952 : BOOST_PP_SUB(DIM, 1), \
953 : (EQUATION_OF_STATE_FUNCTIONS_1D, EQUATION_OF_STATE_FUNCTIONS_2D, \
954 : EQUATION_OF_STATE_FUNCTIONS_3D))))
955 :
956 : /// \cond
957 : #define EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS_HELPER(r, DIM, \
958 : FUNCTION_NAME) \
959 : template <class DataType> \
960 : Scalar<DataType> FUNCTION_NAME##_impl(BOOST_PP_REPEAT( \
961 : DIM, EQUATION_OF_STATE_ARGUMENTS_EXPAND, DataType)) const;
962 : /// \endcond
963 :
964 0 : #define EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(DIM) \
965 : BOOST_PP_LIST_FOR_EACH( \
966 : EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS_HELPER, DIM, \
967 : BOOST_PP_TUPLE_TO_LIST(BOOST_PP_TUPLE_ELEM( \
968 : BOOST_PP_SUB(DIM, 1), \
969 : (EQUATION_OF_STATE_FUNCTIONS_1D, EQUATION_OF_STATE_FUNCTIONS_2D, \
970 : EQUATION_OF_STATE_FUNCTIONS_3D))))
|
