Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <pup.h>
7 :
8 : #include "DataStructures/Tensor/TypeAliases.hpp"
9 : #include "Options/String.hpp"
10 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
11 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
12 : #include "PointwiseFunctions/AnalyticSolutions/GrMhd/Solutions.hpp"
13 : #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp"
14 : #include "PointwiseFunctions/Hydro/Tags.hpp"
15 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
16 : #include "Utilities/Requires.hpp"
17 : #include "Utilities/Serialization/CharmPupable.hpp"
18 : #include "Utilities/TMPL.hpp"
19 : #include "Utilities/TaggedTuple.hpp"
20 :
21 : namespace grmhd::Solutions {
22 :
23 : /*!
24 : * \brief Bondi-Michel accretion \cite Michel1972 with superposed magnetic field
25 : * in Schwarzschild spacetime in Cartesian Kerr-Schild coordinates.
26 : *
27 : * An analytic solution to the 3-D GRMHD system. The user specifies the sonic
28 : * radius \f$r_c\f$ and sonic rest mass density \f$\rho_c\f$, which are the
29 : * radius and rest mass density at the sonic point, the radius at which the
30 : * fluid's Eulerian velocity as seen by a distant observer overtakes the local
31 : * sound speed \f$c_{s,c}\f$. With a specified polytropic exponent \f$\gamma\f$,
32 : * these quantities can be related to the sound speed at infinity
33 : * \f$c_{s,\infty}\f$ using the following relations:
34 : *
35 : * \f{align*}
36 : * c_{s,c}^2 &= \frac{M}{2r_c - 3M} \\
37 : * c_{s,\infty}^2 &= \gamma - 1 + (c_{s,c}^2 - \gamma + 1)\sqrt{1 + 3c_{s,c}^2}
38 : * \f}
39 : *
40 : * In the case of the interstellar medium, the sound
41 : * speed is \f$\approx 10^{-4}\f$, which results in a sonic radius of
42 : * \f$\approx 10^8 M\f$ for \f$\gamma \neq 5/3\f$ \cite RezzollaBook.
43 : *
44 : * The density is found via root-finding, through the
45 : * Bernoulli equation. As one approaches the sonic radius, a second root makes
46 : * an appearance and one must take care to bracket the correct root. This is
47 : * done by using the upper bound
48 : * \f$\frac{\dot{M}}{4\pi}\sqrt{\frac{2}{Mr^3}}\f$.
49 : *
50 : * Additionally specified by the user are the polytropic exponent \f$\gamma\f$,
51 : * and the strength parameter of the magnetic field \f$B_0\f$.
52 : * In Cartesian Kerr-Schild coordinates \f$(x, y, z)\f$, where
53 : * \f$ r = \sqrt{x^2 + y^2 + z^2}\f$, the superposed magnetic field is
54 : * \cite Etienne2010ui
55 : *
56 : * \f{align*}
57 : * B^i(\vec{x},t) = \frac{B_0 M^2}{r^3 \sqrt{\gamma}}x^i
58 : * =\frac{B_0 M^2}{r^3 \sqrt{1 + 2M/r}}x^i.
59 : * \f}
60 : *
61 : * The accretion rate is
62 : *
63 : * \f{align*}{
64 : * \dot{M}=4\pi r^2\rho u^r,
65 : * \f}
66 : *
67 : * and at the sonic radius
68 : *
69 : * \f{align*}{
70 : * \dot{M}_c=\sqrt{8}\pi \sqrt{M}r_c^{3/2}\rho_c.
71 : * \f}
72 : *
73 : * The polytropic constant is given by
74 : *
75 : * \f{align*}{
76 : * K=\frac{1}{\gamma\rho_c^{\gamma-1}}
77 : * \left[\frac{M(\gamma-1)}{(2r_c-3M)(\gamma-1)-M}\right].
78 : * \f}
79 : *
80 : * The density as a function of the sound speed is
81 : *
82 : * \f{align*}{
83 : * \rho^{\gamma-1}=\frac{(\gamma-1)c_s^2}{\gamma K(\gamma-1-c_s^2)}.
84 : * \f}
85 : *
86 : * ## Horizon quantities, \f$\gamma\ne5/3\f$
87 : *
88 : * The density at the horizon is given by:
89 : *
90 : * \f{align*}{
91 : * \rho_h\simeq\frac{1}{16}
92 : * \left(\frac{5-3\gamma}{2}\right)^{(3\gamma-5)/[2(\gamma-1)]}
93 : * \frac{\rho_\infty}{c_{s,\infty}^3}.
94 : * \f}
95 : *
96 : * Using the Lorentz invariance of \f$b^2\f$ we evaluate:
97 : *
98 : * \f{align*}{
99 : * b^2=\frac{B^2}{W^2}+(B^i v_i)^2=
100 : * B^r B^r(1-\gamma_{rr}v^r v^r)+B^r B^r v^r v^r
101 : * =B^r B^r = \frac{B_0^2 M^4}{r^4},
102 : * \f}
103 : *
104 : * where \f$r\f$ is the Cartesian Kerr-Schild radius, which is equal to the
105 : * areal radius for a non-spinning black hole. At the horizon we get
106 : *
107 : * \f{align*}{
108 : * b^2_h=\frac{B^2_0}{16}.
109 : * \f}
110 : *
111 : * Finally, we get
112 : *
113 : * \f{align*}{
114 : * B_0 = 4 \sqrt{b^2_h} = 4\sqrt{\rho_h} \sqrt{\frac{b^2_h}{\rho_h}},
115 : * \f}
116 : *
117 : * where the last equality is useful for comparison to papers that give
118 : * \f$b^2_h/\rho_h\f$.
119 : *
120 : * To help with comparing to other codes the following script can be used to
121 : * compute \f$b^2_h/\rho_h\f$:
122 : *
123 : * \code{.py}
124 : * #!/bin/env python
125 : *
126 : * import numpy as np
127 : *
128 : * # Input parameters
129 : * B_0 = 18
130 : * r_c = 8
131 : * rho_c = 1 / 16
132 : * gamma = 4 / 3
133 : * mass = 1
134 : *
135 : * K = 1 / (gamma * rho_c**(gamma - 1)) * ((gamma - 1) * mass) / (
136 : * (2 * r_c - 3 * mass) * (gamma - 1) - mass)
137 : * c_s_c = mass / (2 * r_c - 3 * mass)
138 : * c_inf = gamma - 1 + (c_s_c - gamma + 1) * np.sqrt(1. + 3. * c_s_c)
139 : * rho_inf = ((gamma - 1) * c_inf / (gamma * K *
140 : * (gamma - 1 - c_inf)))**(1. / (gamma - 1.))
141 : * rho_h = 1. / 16. * (2.5 - 1.5 * gamma)**(
142 : * (3 * gamma - 5) / (2 * (gamma - 1))) * rho_inf / (c_inf**1.5)
143 : *
144 : * print("B_0", B_0)
145 : * print("r_c: ", r_c)
146 : * print("rho_c", rho_c)
147 : * print("b_h^2/rho_h: ", B_0**2 / (16. * rho_h))
148 : * print("gamma: ", gamma)
149 : * \endcode
150 : *
151 : * ## Horizon quantities, \f$\gamma=5/3\f$
152 : *
153 : * The density at the horizon is given by:
154 : *
155 : * \f{align*}{
156 : * \rho_h\simeq \frac{1}{16}\frac{\rho_\infty}{u_h c_{s,\infty}^3},
157 : * \f}
158 : *
159 : * which gives \cite RezzollaBook
160 : *
161 : * \f{align*}{
162 : * \rho_h\simeq 0.08\frac{\rho_\infty}{c_{s,\infty}^3}.
163 : * \f}
164 : *
165 : * The magnetic field \f$b^2\f$ is the same as the \f$\gamma\ne5/3\f$.
166 : */
167 1 : class BondiMichel : public virtual evolution::initial_data::InitialData,
168 : public AnalyticSolution,
169 : public MarkAsAnalyticSolution {
170 : protected:
171 : template <typename DataType>
172 : struct IntermediateVars;
173 :
174 : public:
175 0 : using equation_of_state_type = EquationsOfState::PolytropicFluid<true>;
176 :
177 : /// The mass of the black hole.
178 1 : struct Mass {
179 0 : using type = double;
180 0 : static constexpr Options::String help = {"Mass of the black hole."};
181 0 : static type lower_bound() { return 0.0; }
182 : };
183 :
184 : /// The radius at which the fluid becomes supersonic.
185 1 : struct SonicRadius {
186 0 : using type = double;
187 0 : static constexpr Options::String help = {
188 : "Schwarzschild radius where fluid becomes supersonic."};
189 0 : static type lower_bound() { return 0.0; }
190 : };
191 :
192 : /// The rest mass density of the fluid at the sonic radius.
193 1 : struct SonicDensity {
194 0 : using type = double;
195 0 : static constexpr Options::String help = {
196 : "The density of the fluid at the sonic radius."};
197 0 : static type lower_bound() { return 0.0; }
198 : };
199 :
200 : /// The polytropic exponent for the polytropic fluid.
201 1 : struct PolytropicExponent {
202 0 : using type = double;
203 0 : static constexpr Options::String help = {
204 : "The polytropic exponent for the polytropic fluid."};
205 0 : static type lower_bound() { return 1.0; }
206 : };
207 :
208 : /// The strength of the radial magnetic field.
209 1 : struct MagFieldStrength {
210 0 : using type = double;
211 0 : static constexpr Options::String help = {
212 : "The strength of the radial magnetic field."};
213 : };
214 :
215 0 : using options = tmpl::list<Mass, SonicRadius, SonicDensity,
216 : PolytropicExponent, MagFieldStrength>;
217 0 : static constexpr Options::String help = {
218 : "Bondi-Michel solution with a radial magnetic field using \n"
219 : "the Schwarzschild coordinate system. Quantities prefixed with \n"
220 : "`sonic` refer to field quantities evaluated at the radius \n"
221 : "where the fluid speed overtakes the sound speed."};
222 :
223 0 : BondiMichel() = default;
224 0 : BondiMichel(const BondiMichel& /*rhs*/) = default;
225 0 : BondiMichel& operator=(const BondiMichel& /*rhs*/) = default;
226 0 : BondiMichel(BondiMichel&& /*rhs*/) = default;
227 0 : BondiMichel& operator=(BondiMichel&& /*rhs*/) = default;
228 0 : ~BondiMichel() override = default;
229 :
230 0 : BondiMichel(double mass, double sonic_radius, double sonic_density,
231 : double polytropic_exponent, double mag_field_strength);
232 :
233 0 : auto get_clone() const
234 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
235 :
236 : /// \cond
237 : explicit BondiMichel(CkMigrateMessage* msg);
238 : using PUP::able::register_constructor;
239 : WRAPPED_PUPable_decl_template(BondiMichel);
240 : /// \endcond
241 :
242 : /// Retrieve a collection of hydro variables at `(x, t)`
243 : template <typename DataType, typename... Tags>
244 1 : tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
245 : const double /*t*/,
246 : tmpl::list<Tags...> /*meta*/) const {
247 : // non-const so we can move out the metric vars. We assume that no variable
248 : // is being retrieved more than once, which would cause problems with
249 : // TaggedTuple anyway.
250 : auto intermediate_vars = IntermediateVars<DataType>{
251 : rest_mass_density_at_infinity_,
252 : mass_accretion_rate_over_four_pi_,
253 : mass_,
254 : polytropic_constant_,
255 : polytropic_exponent_,
256 : bernoulli_constant_squared_minus_one_,
257 : sonic_radius_,
258 : sonic_density_,
259 : x,
260 : tmpl2::flat_any_v<
261 : not tmpl::list_contains_v<hydro::grmhd_tags<DataType>, Tags>...>,
262 : background_spacetime_};
263 : return {get<Tags>(variables(x, tmpl::list<Tags>{}, intermediate_vars))...};
264 : }
265 :
266 : template <typename DataType, typename Tag>
267 0 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
268 : const double /*t*/, // NOLINT
269 : tmpl::list<Tag> /*meta*/) const {
270 : return variables(
271 : x, tmpl::list<Tag>{},
272 : IntermediateVars<DataType>{
273 : rest_mass_density_at_infinity_, mass_accretion_rate_over_four_pi_,
274 : mass_, polytropic_constant_, polytropic_exponent_,
275 : bernoulli_constant_squared_minus_one_, sonic_radius_,
276 : sonic_density_, x,
277 : not tmpl::list_contains_v<hydro::grmhd_tags<DataType>, Tag>,
278 : background_spacetime_});
279 : }
280 :
281 : // NOLINTNEXTLINE(google-runtime-references)
282 0 : void pup(PUP::er& /*p*/) override;
283 0 : const EquationsOfState::PolytropicFluid<true>& equation_of_state() const {
284 : return equation_of_state_;
285 : }
286 :
287 : protected:
288 0 : friend bool operator==(const BondiMichel& lhs, const BondiMichel& rhs);
289 :
290 : template <typename DataType>
291 0 : auto variables(const tnsr::I<DataType, 3>& x,
292 : tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/,
293 : const IntermediateVars<DataType>& vars) const
294 : -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
295 :
296 : template <typename DataType>
297 0 : auto variables(const tnsr::I<DataType, 3>& x,
298 : tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/,
299 : const IntermediateVars<DataType>& vars) const
300 : -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
301 :
302 : template <typename DataType>
303 0 : auto variables(
304 : const tnsr::I<DataType, 3>& x,
305 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/,
306 : const IntermediateVars<DataType>& vars) const
307 : -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
308 :
309 : template <typename DataType>
310 0 : auto variables(const tnsr::I<DataType, 3>& x,
311 : tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/,
312 : const IntermediateVars<DataType>& vars) const
313 : -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
314 :
315 : template <typename DataType>
316 0 : auto variables(const tnsr::I<DataType, 3>& x,
317 : tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/,
318 : const IntermediateVars<DataType>& vars) const
319 : -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>>;
320 :
321 : template <typename DataType>
322 0 : auto variables(const tnsr::I<DataType, 3>& x,
323 : tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/,
324 : const IntermediateVars<DataType>& vars) const
325 : -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
326 :
327 : template <typename DataType>
328 0 : auto variables(const tnsr::I<DataType, 3>& x,
329 : tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/,
330 : const IntermediateVars<DataType>& vars) const
331 : -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
332 :
333 : template <typename DataType>
334 0 : auto variables(
335 : const tnsr::I<DataType, 3>& x,
336 : tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/,
337 : const IntermediateVars<DataType>& vars) const
338 : -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
339 :
340 : template <typename DataType>
341 0 : auto variables(const tnsr::I<DataType, 3>& x,
342 : tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/,
343 : const IntermediateVars<DataType>& vars) const
344 : -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
345 :
346 : template <typename DataType>
347 0 : auto variables(const tnsr::I<DataType, 3>& x,
348 : tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/,
349 : const IntermediateVars<DataType>& vars) const
350 : -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
351 :
352 : template <typename DataType, typename Tag,
353 : Requires<not tmpl::list_contains_v<
354 : tmpl::push_back<hydro::grmhd_tags<DataType>,
355 : hydro::Tags::SpecificEnthalpy<DataType>>,
356 : Tag>> = nullptr>
357 0 : tuples::TaggedTuple<Tag> variables(
358 : const tnsr::I<DataType, 3>& /*x*/, tmpl::list<Tag> /*meta*/,
359 : const IntermediateVars<DataType>& vars) const {
360 : return {std::move(get<Tag>(vars.kerr_schild_soln))};
361 : }
362 :
363 : template <typename DataType>
364 0 : struct IntermediateVars {
365 0 : IntermediateVars(double rest_mass_density_at_infinity,
366 : double in_mass_accretion_rate_over_four_pi, double in_mass,
367 : double in_polytropic_constant,
368 : double in_polytropic_exponent,
369 : double in_bernoulli_constant_squared_minus_one,
370 : double in_sonic_radius, double in_sonic_density,
371 : const tnsr::I<DataType, 3>& x, bool need_spacetime,
372 : const gr::Solutions::KerrSchild& background_spacetime);
373 0 : DataType radius{};
374 0 : DataType rest_mass_density{};
375 0 : DataType electron_fraction{};
376 0 : double mass_accretion_rate_over_four_pi{};
377 0 : double mass{};
378 0 : double polytropic_constant{};
379 0 : double polytropic_exponent{};
380 0 : double bernoulli_constant_squared_minus_one{};
381 0 : double sonic_radius{};
382 0 : double sonic_density{};
383 0 : double bernoulli_root_function(double rest_mass_density_guess,
384 : double current_radius) const;
385 : tuples::tagged_tuple_from_typelist<
386 : typename gr::Solutions::KerrSchild::tags<DataType>>
387 0 : kerr_schild_soln{};
388 : };
389 :
390 0 : double mass_ = std::numeric_limits<double>::signaling_NaN();
391 0 : double sonic_radius_ = std::numeric_limits<double>::signaling_NaN();
392 0 : double sonic_density_ = std::numeric_limits<double>::signaling_NaN();
393 0 : double polytropic_exponent_ = std::numeric_limits<double>::signaling_NaN();
394 0 : double mag_field_strength_ = std::numeric_limits<double>::signaling_NaN();
395 0 : double sonic_fluid_speed_squared_ =
396 : std::numeric_limits<double>::signaling_NaN();
397 0 : double sonic_sound_speed_squared_ =
398 : std::numeric_limits<double>::signaling_NaN();
399 0 : double polytropic_constant_ = std::numeric_limits<double>::signaling_NaN();
400 0 : double mass_accretion_rate_over_four_pi_ =
401 : std::numeric_limits<double>::signaling_NaN();
402 0 : double bernoulli_constant_squared_minus_one_ =
403 : std::numeric_limits<double>::signaling_NaN();
404 0 : double rest_mass_density_at_infinity_ =
405 : std::numeric_limits<double>::signaling_NaN();
406 0 : EquationsOfState::PolytropicFluid<true> equation_of_state_{};
407 0 : gr::Solutions::KerrSchild background_spacetime_{};
408 : };
409 :
410 0 : bool operator!=(const BondiMichel& lhs, const BondiMichel& rhs);
411 :
412 : } // namespace grmhd::Solutions
|
