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/AnalyticSolutions/AnalyticSolution.hpp"
12 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp"
13 : #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Solutions.hpp"
14 : #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp" // IWYU pragma: keep
15 : #include "PointwiseFunctions/Hydro/Tags.hpp"
16 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
17 : #include "PointwiseFunctions/Hydro/Temperature.hpp"
18 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
19 : #include "Utilities/ForceInline.hpp"
20 : #include "Utilities/Requires.hpp"
21 : #include "Utilities/Serialization/CharmPupable.hpp"
22 : #include "Utilities/TMPL.hpp"
23 : #include "Utilities/TaggedTuple.hpp"
24 :
25 : /// \cond
26 : namespace PUP {
27 : class er; // IWYU pragma: keep
28 : } // namespace PUP
29 : namespace grmhd::AnalyticData {
30 : class MagnetizedFmDisk;
31 : } // namespace grmhd::AnalyticData
32 : /// \endcond
33 :
34 : // IWYU pragma: no_include <pup.h>
35 :
36 : namespace RelativisticEuler {
37 1 : namespace Solutions {
38 :
39 : /*!
40 : * \brief Fluid disk orbiting a Kerr black hole
41 : *
42 : * The Fishbone-Moncrief solution to the 3D relativistic Euler system
43 : * \cite Fishbone1976apj, representing the isentropic flow of a thick fluid disk
44 : * orbiting a Kerr black hole. In Boyer-Lindquist coordinates
45 : * \f$(t, r, \theta, \phi)\f$, the flow is assumed to be purely toroidal,
46 : *
47 : * \f{align*}
48 : * u^\mu = (u^t, 0, 0, u^\phi),
49 : * \f}
50 : *
51 : * where \f$u^\mu\f$ is the 4-velocity. Then, all the fluid quantities are
52 : * assumed to share the same symmetries as those of the background spacetime,
53 : * namely they are stationary (independent of \f$t\f$), and axially symmetric
54 : * (independent of \f$\phi\f$).
55 : *
56 : * Self-gravity is neglected, so that the fluid
57 : * variables are determined as functions of the metric. Following the treatment
58 : * by Kozlowski et al. \cite Kozlowski1978aa (but using signature +2)
59 : * the solution is expressed in terms of the quantities
60 : *
61 : * \f{align*}
62 : * \Omega &= \dfrac{u^\phi}{u^t},\\
63 : * W &= W_\text{in} - \int_{p_\text{in}}^p\frac{dp}{e + p},
64 : * \f}
65 : *
66 : * where \f$\Omega\f$ is the angular velocity, \f$p\f$ is the fluid pressure,
67 : * \f$e\f$ is the energy density, and \f$W\f$ is an auxiliary quantity
68 : * interpreted in the Newtonian limit as the total (gravitational + centrifugal)
69 : * potential. \f$W_\text{in}\f$ and \f$p_\text{in}\f$ are the potential and
70 : * the pressure at the radius of the inner edge, i.e. the closest edge to the
71 : * black hole. Here we assume \f$p_\text{in} = 0.\f$ The solution to the Euler
72 : * equation is then
73 : *
74 : * \f{align*}
75 : * (u^t)^2 &= \frac{A}{2\Delta\Sigma}\left(1 +
76 : * \sqrt{1 + \frac{4l^2\Delta \Sigma^2}{A^2\sin^2\theta}}\right)\\
77 : * \Omega &= \frac{\Sigma}{A (u^t)^2}\frac{l}{\sin^2\theta} + \frac{2Mra}{A}\\
78 : * u^\phi &= \Omega u^t\\
79 : * W &= l\Omega - \ln u^t,
80 : * \f}
81 : *
82 : * where
83 : *
84 : * \f{align*}
85 : * \Sigma = r^2 + a^2\cos^2\theta\qquad
86 : * \Delta = r^2 - 2Mr + a^2\qquad
87 : * A = (r^2 + a^2)^2 - \Delta a^2 \sin^2\theta
88 : * \f}
89 : *
90 : * and \f$l = u_\phi u^t\f$ is the so-called angular momentum per unit
91 : * intertial mass, which is a parameter defining an
92 : * individual disk. In deriving the solution, an integration constant has been
93 : * chosen so that \f$ W\longrightarrow 0\f$ as \f$r\longrightarrow \infty\f$,
94 : * in accordance with the Newtonian limit. Note that, from its definition,
95 : * equipotential contours coincide with isobaric contours. Physically, the
96 : * matter can fill each of the closed surfaces \f$W = \text{const}\f$, giving
97 : * rise to an orbiting thick disk. For \f$W > 0\f$, all equipotentials are open,
98 : * whereas for \f$W < 0\f$, some of them will be closed. Should a disk exist,
99 : * the pressure reaches a maximum value on the equator at a coordinate radius
100 : * \f$r_\text{max}\f$ that is related to the angular momentum per unit inertial
101 : * mass via
102 : *
103 : * \f{align*}
104 : * l = \dfrac{M^{1/2}(r_\text{max}^{3/2} + aM^{1/2})(a^2 - 2aM^{1/2}
105 : * r_\text{max}^{1/2} + r_\text{max}^2)}{2aM^{1/2}r_\text{max}^{3/2} +
106 : * (r_\text{max} - 3M)r_\text{max}^2}.
107 : * \f}
108 : *
109 : * Once \f$W\f$ is determined, an equation of state is required in order to
110 : * obtain the thermodynamic variables. If the flow is isentropic, the specific
111 : * enthalpy can readily be obtained from the first and second laws of
112 : * thermodynamics: one has
113 : *
114 : * \f{align*}
115 : * \frac{dp}{e + p} = \frac{dh}{h}
116 : * \f}
117 : *
118 : * so that
119 : *
120 : * \f{align*}
121 : * h = h_\text{in}\exp(W_\text{in} - W),
122 : * \f}
123 : *
124 : * and the pressure can be obtained from a thermodynamic relation of the form
125 : * \f$h = h(p)\f$. Here we assume a polytropic relation
126 : *
127 : * \f{align*}
128 : * p = K\rho^\gamma.
129 : * \f}
130 : *
131 : * Once all the variables are known in Boyer-Lindquist (or Kerr) coordinates, it
132 : * is straightforward to write them in Cartesian Kerr-Schild coordinates. The
133 : * coordinate transformation in gr::KerrSchildCoords helps read the Jacobian
134 : * matrix, which, applied to the azimuthal flow of the disk, gives
135 : *
136 : * \f{align*}
137 : * u_\text{KS}^\mu = u^t(1, -y\Omega, x\Omega, 0),
138 : * \f}
139 : *
140 : * where \f$u^t\f$ and \f$\Omega\f$ are now understood as functions of the
141 : * Kerr-Schild coordinates. Finally, the spatial velocity can be readily
142 : * obtained from its definition,
143 : *
144 : * \f{align*}
145 : * \alpha v^i = \frac{u^i}{u^t} + \beta^i,
146 : * \f}
147 : *
148 : * where \f$\alpha\f$ and \f$\beta^i\f$ are the lapse and the shift,
149 : * respectively.
150 : *
151 : * \note Kozlowski et al. \cite Kozlowski1978aa denote
152 : * \f$l_* = u_\phi u^t\f$ in order to
153 : * distinguish this quantity from their own definition \f$l = - u_\phi/u_t\f$.
154 : *
155 : * \note When using Kerr-Schild coordinates, the horizon that is at
156 : * constant \f$r\f$ is not spherical, but instead spheroidal. This could make
157 : * application of boundary condition and computing various fluxes
158 : * across the horizon more complicated than they need to be.
159 : * Thus, we use Spherical Kerr-Schild coordinates,
160 : * see gr::Solutions::SphericalKerrSchild, in which constant \f$r\f$
161 : * is spherical. Because we compute variables in Kerr-Schild coordinates,
162 : * there is a necessary extra step of transforming them back to
163 : * Spherical Kerr-Schild coordinates.
164 : *
165 : */
166 1 : class FishboneMoncriefDisk
167 : : public virtual evolution::initial_data::InitialData,
168 : public MarkAsAnalyticSolution,
169 : public AnalyticSolution<3>,
170 : public hydro::TemperatureInitialization<FishboneMoncriefDisk> {
171 : protected:
172 : template <typename DataType>
173 : struct IntermediateVariables;
174 :
175 : public:
176 0 : using equation_of_state_type = EquationsOfState::PolytropicFluid<true>;
177 :
178 : /// The mass of the black hole, \f$M\f$.
179 1 : struct BhMass {
180 0 : using type = double;
181 0 : static constexpr Options::String help = {"The mass of the black hole."};
182 0 : static type lower_bound() { return 0.0; }
183 : };
184 : /// The dimensionless black hole spin, \f$\chi = a/M\f$.
185 1 : struct BhDimlessSpin {
186 0 : using type = double;
187 0 : static constexpr Options::String help = {
188 : "The dimensionless black hole spin."};
189 0 : static type lower_bound() { return 0.0; }
190 0 : static type upper_bound() { return 1.0; }
191 : };
192 : /// The radial coordinate of the inner edge of the disk, in units of \f$M\f$.
193 1 : struct InnerEdgeRadius {
194 0 : using type = double;
195 0 : static constexpr Options::String help = {
196 : "The radial coordinate of the inner edge of the disk."};
197 : };
198 : /// The radial coordinate of the maximum pressure, in units of \f$M\f$.
199 1 : struct MaxPressureRadius {
200 0 : using type = double;
201 0 : static constexpr Options::String help = {
202 : "The radial coordinate of the maximum pressure."};
203 : };
204 : /// The polytropic constant of the fluid.
205 1 : struct PolytropicConstant {
206 0 : using type = double;
207 0 : static constexpr Options::String help = {
208 : "The polytropic constant of the fluid."};
209 0 : static type lower_bound() { return 0.; }
210 : };
211 : /// The polytropic exponent of the fluid.
212 1 : struct PolytropicExponent {
213 0 : using type = double;
214 0 : static constexpr Options::String help = {
215 : "The polytropic exponent of the fluid."};
216 0 : static type lower_bound() { return 1.; }
217 : };
218 :
219 0 : using options =
220 : tmpl::list<BhMass, BhDimlessSpin, InnerEdgeRadius, MaxPressureRadius,
221 : PolytropicConstant, PolytropicExponent>;
222 0 : static constexpr Options::String help = {
223 : "Fluid disk orbiting a Kerr black hole."};
224 :
225 0 : FishboneMoncriefDisk() = default;
226 0 : FishboneMoncriefDisk(const FishboneMoncriefDisk& /*rhs*/) = default;
227 0 : FishboneMoncriefDisk& operator=(const FishboneMoncriefDisk& /*rhs*/) =
228 : default;
229 0 : FishboneMoncriefDisk(FishboneMoncriefDisk&& /*rhs*/) = default;
230 0 : FishboneMoncriefDisk& operator=(FishboneMoncriefDisk&& /*rhs*/) = default;
231 0 : ~FishboneMoncriefDisk() override = default;
232 :
233 0 : FishboneMoncriefDisk(double bh_mass, double bh_dimless_spin,
234 : double inner_edge_radius, double max_pressure_radius,
235 : double polytropic_constant, double polytropic_exponent);
236 :
237 0 : auto get_clone() const
238 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
239 :
240 : /// \cond
241 : explicit FishboneMoncriefDisk(CkMigrateMessage* msg);
242 : using PUP::able::register_constructor;
243 : WRAPPED_PUPable_decl_template(FishboneMoncriefDisk);
244 : /// \endcond
245 :
246 : // Eventually, if we implement a gr::Solutions::BoyerLindquist
247 : // black hole, the following two functions aren't needed, and the algebra
248 : // of the three functions after can be simplified by using the corresponding
249 : // lapse, shift, and spatial metric.
250 : template <typename DataType>
251 0 : DataType sigma(const DataType& r_sqrd, const DataType& sin_theta_sqrd) const;
252 :
253 : template <typename DataType>
254 0 : DataType inv_ucase_a(const DataType& r_sqrd, const DataType& sin_theta_sqrd,
255 : const DataType& delta) const;
256 :
257 : template <typename DataType>
258 0 : DataType four_velocity_t_sqrd(const DataType& r_sqrd,
259 : const DataType& sin_theta_sqrd) const;
260 :
261 : template <typename DataType>
262 0 : DataType angular_velocity(const DataType& r_sqrd,
263 : const DataType& sin_theta_sqrd) const;
264 :
265 : template <typename DataType>
266 0 : DataType potential(const DataType& r_sqrd,
267 : const DataType& sin_theta_sqrd) const;
268 :
269 : template <typename DataType>
270 0 : using tags =
271 : tmpl::append<hydro::grmhd_tags<DataType>,
272 : typename gr::Solutions::SphericalKerrSchild::tags<DataType>>;
273 :
274 : /// @{
275 : /// The variables in Cartesian Spherical-Kerr-Schild coordinates at `(x, t)`
276 : template <typename DataType, typename... Tags>
277 1 : tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
278 : const double /*t*/,
279 : tmpl::list<Tags...> /*meta*/) const {
280 : // Can't store IntermediateVariables as member variable because we need to
281 : // be threadsafe.
282 : IntermediateVariables<DataType> vars(x);
283 : return {std::move(
284 : get<Tags>(variables(x, tmpl::list<Tags>{}, make_not_null(&vars))))...};
285 : }
286 :
287 : template <typename DataType, typename Tag>
288 1 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
289 : const double /*t*/,
290 : tmpl::list<Tag> /*meta*/) const {
291 : // Can't store IntermediateVariables as member variable because we need to
292 : // be threadsafe.
293 : IntermediateVariables<DataType> vars(x);
294 : return variables(x, tmpl::list<Tag>{}, make_not_null(&vars));
295 : }
296 : /// @}
297 :
298 : // NOLINTNEXTLINE(google-runtime-references)
299 0 : void pup(PUP::er& p) override;
300 :
301 0 : const EquationsOfState::PolytropicFluid<true>& equation_of_state() const {
302 : return equation_of_state_;
303 : }
304 :
305 : protected:
306 : template <typename DataType>
307 0 : auto variables(const tnsr::I<DataType, 3>& x,
308 : tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/,
309 : gsl::not_null<IntermediateVariables<DataType>*> vars) const
310 : -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
311 :
312 : template <typename DataType>
313 0 : auto variables(const tnsr::I<DataType, 3>& x,
314 : tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/,
315 : gsl::not_null<IntermediateVariables<DataType>*> vars) const
316 : -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
317 :
318 : template <typename DataType>
319 0 : auto variables(const tnsr::I<DataType, 3>& x,
320 : tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/,
321 : gsl::not_null<IntermediateVariables<DataType>*> vars) const
322 : -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
323 :
324 : template <typename DataType>
325 0 : auto variables(const tnsr::I<DataType, 3>& x,
326 : tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/,
327 : gsl::not_null<IntermediateVariables<DataType>*> vars) const
328 : -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
329 :
330 : template <typename DataType>
331 0 : auto variables(const tnsr::I<DataType, 3>& x,
332 : tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/,
333 : gsl::not_null<IntermediateVariables<DataType>*> vars) const
334 : -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>>;
335 :
336 : template <typename DataType>
337 0 : auto variables(
338 : const tnsr::I<DataType, 3>& x,
339 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/,
340 : gsl::not_null<IntermediateVariables<DataType>*> vars) const
341 : -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
342 :
343 : template <typename DataType>
344 0 : auto variables(const tnsr::I<DataType, 3>& x,
345 : tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/,
346 : gsl::not_null<IntermediateVariables<DataType>*> vars) const
347 : -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
348 :
349 : template <typename DataType>
350 0 : auto variables(const tnsr::I<DataType, 3>& x,
351 : tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/,
352 : gsl::not_null<IntermediateVariables<DataType>*> vars) const
353 : -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
354 :
355 : template <typename DataType>
356 0 : auto variables(const tnsr::I<DataType, 3>& x,
357 : tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/,
358 : gsl::not_null<IntermediateVariables<DataType>*> vars) const
359 : -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
360 :
361 : template <typename DataType>
362 0 : auto variables(
363 : const tnsr::I<DataType, 3>& x,
364 : tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/,
365 : gsl::not_null<IntermediateVariables<DataType>*> vars) const
366 : -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
367 :
368 : // Grab the metric variables
369 : template <typename DataType, typename Tag,
370 : Requires<not tmpl::list_contains_v<
371 : tmpl::push_back<hydro::grmhd_tags<DataType>,
372 : hydro::Tags::SpecificEnthalpy<DataType>,
373 : hydro::Tags::SpatialVelocity<DataType, 3>,
374 : hydro::Tags::LorentzFactor<DataType>>,
375 : Tag>> = nullptr>
376 0 : tuples::TaggedTuple<Tag> variables(
377 : const tnsr::I<DataType, 3>& x, tmpl::list<Tag> /*meta*/,
378 : gsl::not_null<IntermediateVariables<DataType>*> vars) const {
379 : return {get<Tag>(background_spacetime_.variables(
380 : x, 0.0, tmpl::list<Tag>{},
381 : make_not_null(&vars->sph_kerr_schild_cache)))};
382 : }
383 :
384 : template <typename DataType, typename Func>
385 0 : void variables_impl(gsl::not_null<IntermediateVariables<DataType>*> vars,
386 : Func f) const;
387 :
388 : // Intermediate variables needed to set several of the Fishbone-Moncrief
389 : // solution's variables.
390 :
391 : template <typename DataType>
392 0 : struct IntermediateVariables {
393 0 : explicit IntermediateVariables(const tnsr::I<DataType, 3>& x);
394 :
395 0 : DataType r_squared{};
396 0 : DataType sin_theta_squared{};
397 : gr::Solutions::SphericalKerrSchild::IntermediateVars<DataType,
398 : Frame::Inertial>
399 0 : sph_kerr_schild_cache =
400 : gr::Solutions::SphericalKerrSchild::IntermediateVars<
401 : DataType, Frame::Inertial>(0);
402 : };
403 :
404 0 : friend bool operator==(const FishboneMoncriefDisk& lhs,
405 : const FishboneMoncriefDisk& rhs);
406 0 : friend class grmhd::AnalyticData::MagnetizedFmDisk;
407 :
408 0 : double bh_mass_ = std::numeric_limits<double>::signaling_NaN();
409 0 : double bh_spin_a_ = std::numeric_limits<double>::signaling_NaN();
410 0 : double inner_edge_radius_ = std::numeric_limits<double>::signaling_NaN();
411 0 : double max_pressure_radius_ = std::numeric_limits<double>::signaling_NaN();
412 0 : double polytropic_constant_ = std::numeric_limits<double>::signaling_NaN();
413 0 : double polytropic_exponent_ = std::numeric_limits<double>::signaling_NaN();
414 0 : double angular_momentum_ = std::numeric_limits<double>::signaling_NaN();
415 0 : EquationsOfState::PolytropicFluid<true> equation_of_state_{};
416 0 : gr::Solutions::SphericalKerrSchild background_spacetime_{};
417 : };
418 :
419 0 : bool operator!=(const FishboneMoncriefDisk& lhs,
420 : const FishboneMoncriefDisk& rhs);
421 : } // namespace Solutions
422 : } // namespace RelativisticEuler
|
