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