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