Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <limits>
8 :
9 : #include "DataStructures/Tensor/TypeAliases.hpp"
10 : #include "Options/String.hpp"
11 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
12 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
13 : #include "Utilities/GenerateInstantiations.hpp"
14 : #include "Utilities/Gsl.hpp"
15 : #include "Utilities/TMPL.hpp"
16 :
17 : /// \cond
18 : namespace PUP {
19 : class er;
20 : } // namespace PUP
21 : class DataVector;
22 : namespace domain {
23 : namespace Tags {
24 : template <size_t Dim, typename Frame>
25 : struct Coordinates;
26 : } // namespace Tags
27 : } // namespace domain
28 : namespace hydro {
29 : namespace Tags {
30 : template <typename DataType>
31 : struct Pressure;
32 : template <typename DataType>
33 : struct RestMassDensity;
34 : template <typename DataType>
35 : struct SpecificInternalEnergy;
36 : template <typename DataType>
37 : struct Temperature;
38 : template <typename DataType>
39 : struct ElectronFraction;
40 : } // namespace Tags
41 : } // namespace hydro
42 : /// \endcond
43 :
44 : /// Contains all variable fixers.
45 : namespace VariableFixing {
46 :
47 : /// \ingroup VariableFixingGroup
48 : /// \brief Applies a pressure and density floor dependent on the distance
49 : /// to the origin.
50 : ///
51 : /// Applies the floors:
52 : /// \f$\rho(r) \geq \rho_{\mathrm{fl}}(r) = C_\rho r^{k_\rho}\f$
53 : /// and \f$P(r) \geq P_{\mathrm{fl}}(r) = C_p r^{k_p}\f$
54 : /// when \f$ r > r_{min}\f$, where \f$C_\rho\f$ is given by the option
55 : /// `ScaleDensityFloor`, \f$k_\rho\f$ is given by the option
56 : /// `PowerDensityFloor`, \f$C_p\f$ is given by the option
57 : /// `ScalePressureFloor`, \f$k_p\f$ is given by the option
58 : /// `PowerPressureFloor`, and \f$r_{min}\f$ is given by the option
59 : /// `MinimumRadius`.
60 : ///
61 : /// \note In \cite Porth2016rfi, the following floors are applied:
62 : /// \f$\rho(r) \geq \rho_{\mathrm{fl}}(r) = 10^{-5}r^{-3/2}\f$
63 : /// and \f$P(r) \geq P_{\mathrm{fl}}(r) = \frac{1}{3} \times 10^{-7}r^{-5/2}\f$
64 : ///
65 : /// \note For 3d equations of state we treat the product \f$\rho T\f$ as the
66 : /// "pressure" that we floor and then compute the pressure from the
67 : /// temperature. The floor is applied after flooring the density.
68 : template <size_t Dim>
69 1 : class RadiallyFallingFloor {
70 : public:
71 : /// \brief The minimum radius at which to begin applying the floors on the
72 : /// density and pressure.
73 1 : struct MinimumRadius {
74 0 : static constexpr Options::String help =
75 : "The radius at which to begin applying the lower bound.";
76 0 : using type = double;
77 0 : static double lower_bound() { return 0.0; }
78 : };
79 :
80 : /// \brief The scale of the floor of the rest mass density.
81 1 : struct ScaleDensityFloor {
82 0 : static constexpr Options::String help =
83 : "The rest mass density floor at r = 1.";
84 0 : using type = double;
85 0 : static double lower_bound() { return 0.0; }
86 : };
87 :
88 : /// \brief The power of the radius of the floor of the rest mass density.
89 1 : struct PowerDensityFloor {
90 0 : static constexpr Options::String help =
91 : "Radial power for the floor of the rest mass density.";
92 0 : using type = double;
93 : };
94 :
95 : /// \brief The scale of the floor of the pressure.
96 1 : struct ScalePressureFloor {
97 0 : static constexpr Options::String help = "The pressure floor at r = 1.";
98 0 : using type = double;
99 0 : static double lower_bound() { return 0.0; }
100 : };
101 :
102 : /// \brief The power of the radius of the floor of the pressure.
103 1 : struct PowerPressureFloor {
104 0 : static constexpr Options::String help =
105 : "The radial power for the floor of the pressure.";
106 0 : using type = double;
107 : };
108 :
109 0 : using options =
110 : tmpl::list<MinimumRadius, ScaleDensityFloor, PowerDensityFloor,
111 : ScalePressureFloor, PowerPressureFloor>;
112 0 : static constexpr Options::String help = {
113 : "Applies a pressure and density floor dependent on the distance to the "
114 : "origin."};
115 :
116 0 : RadiallyFallingFloor(double minimum_radius_at_which_to_apply_floor,
117 : double rest_mass_density_scale,
118 : double rest_mass_density_power, double pressure_scale,
119 : double pressure_power);
120 :
121 0 : RadiallyFallingFloor() = default;
122 0 : RadiallyFallingFloor(const RadiallyFallingFloor& /*rhs*/) = default;
123 0 : RadiallyFallingFloor& operator=(const RadiallyFallingFloor& /*rhs*/) =
124 : default;
125 0 : RadiallyFallingFloor(RadiallyFallingFloor&& /*rhs*/) = default;
126 0 : RadiallyFallingFloor& operator=(RadiallyFallingFloor&& /*rhs*/) = default;
127 0 : ~RadiallyFallingFloor() = default;
128 :
129 : // NOLINTNEXTLINE(google-runtime-references)
130 0 : void pup(PUP::er& p);
131 :
132 0 : using return_tags =
133 : tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
134 : hydro::Tags::Pressure<DataVector>,
135 : hydro::Tags::SpecificInternalEnergy<DataVector>,
136 : hydro::Tags::Temperature<DataVector>,
137 : hydro::Tags::ElectronFraction<DataVector>>;
138 0 : using argument_tags =
139 : tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
140 : hydro::Tags::GrmhdEquationOfState>;
141 :
142 : template <size_t ThermodynamicDim>
143 0 : void operator()(
144 : gsl::not_null<Scalar<DataVector>*> density,
145 : gsl::not_null<Scalar<DataVector>*> pressure,
146 : gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
147 : gsl::not_null<Scalar<DataVector>*> temperature,
148 : gsl::not_null<Scalar<DataVector>*> electron_fraction,
149 : const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
150 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
151 : equation_of_state) const;
152 :
153 : private:
154 : template <size_t LocalDim>
155 : // NOLINTNEXTLINE(readability-redundant-declaration)
156 0 : friend bool operator==(const RadiallyFallingFloor<LocalDim>& lhs,
157 : const RadiallyFallingFloor<LocalDim>& rhs);
158 :
159 0 : double minimum_radius_at_which_to_apply_floor_{
160 : std::numeric_limits<double>::signaling_NaN()};
161 0 : double rest_mass_density_scale_{std::numeric_limits<double>::signaling_NaN()};
162 0 : double rest_mass_density_power_{std::numeric_limits<double>::signaling_NaN()};
163 0 : double pressure_scale_{std::numeric_limits<double>::signaling_NaN()};
164 0 : double pressure_power_{std::numeric_limits<double>::signaling_NaN()};
165 : };
166 :
167 : template <size_t Dim>
168 0 : bool operator!=(const RadiallyFallingFloor<Dim>& lhs,
169 : const RadiallyFallingFloor<Dim>& rhs);
170 :
171 : } // namespace VariableFixing
|
