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" // IWYU pragma: keep
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 : * The density at the horizon is given by:
88 : *
89 : * \f{align*}{
90 : * \rho_h\simeq\frac{1}{16}
91 : * \left(\frac{5-3\gamma}{2}\right)^{(3\gamma-5)/[2(\gamma-1)]}
92 : * \frac{\rho_\infty}{c_{s,\infty}^3}.
93 : * \f}
94 : *
95 : * Using the Lorentz invariance of \f$b^2\f$ we evaluate:
96 : *
97 : * \f{align*}{
98 : * b^2=\frac{B^2}{W^2}+(B^i v_i)^2=
99 : * B^r B^r(1-\gamma_{rr}v^r v^r)+B^r B^r v^r v^r
100 : * =B^r B^r = \frac{B_0^2 M^4}{r^4},
101 : * \f}
102 : *
103 : * where \f$r\f$ is the Cartesian Kerr-Schild radius, which is equal to the
104 : * areal radius for a non-spinning black hole. At the horizon we get
105 : *
106 : * \f{align*}{
107 : * b^2_h=\frac{B^2_0}{16}.
108 : * \f}
109 : *
110 : * Finally, we get
111 : *
112 : * \f{align*}{
113 : * B_0 = 4 \sqrt{b^2_h} = 4\sqrt{\rho_h} \sqrt{\frac{b^2_h}{\rho_h}},
114 : * \f}
115 : *
116 : * where the last equality is useful for comparison to papers that give
117 : * \f$b^2_h/\rho_h\f$.
118 : *
119 : * To help with comparing to other codes the following script can be used to
120 : * compute \f$b^2_h/\rho_h\f$:
121 : *
122 : * \code{.py}
123 : * #!/bin/env python
124 : *
125 : * import numpy as np
126 : *
127 : * # Input parameters
128 : * B_0 = 18
129 : * r_c = 8
130 : * rho_c = 1 / 16
131 : * gamma = 4 / 3
132 : * mass = 1
133 : *
134 : * K = 1 / (gamma * rho_c**(gamma - 1)) * ((gamma - 1) * mass) / (
135 : * (2 * r_c - 3 * mass) * (gamma - 1) - mass)
136 : * c_s_c = mass / (2 * r_c - 3 * mass)
137 : * c_inf = gamma - 1 + (c_s_c - gamma + 1) * np.sqrt(1. + 3. * c_s_c)
138 : * rho_inf = ((gamma - 1) * c_inf / (gamma * K *
139 : * (gamma - 1 - c_inf)))**(1. / (gamma - 1.))
140 : * rho_h = 1. / 16. * (2.5 - 1.5 * gamma)**(
141 : * (3 * gamma - 5) / (2 * (gamma - 1))) * rho_inf / (c_inf**1.5)
142 : *
143 : * print("B_0", B_0)
144 : * print("r_c: ", r_c)
145 : * print("rho_c", rho_c)
146 : * print("b_h^2/rho_h: ", B_0**2 / (16. * rho_h))
147 : * print("gamma: ", gamma)
148 : * \endcode
149 : *
150 : * #### Horizon quantities, \f$\gamma=5/3\f$
151 : * The density at the horizon is given by:
152 : *
153 : * \f{align*}{
154 : * \rho_h\simeq \frac{1}{16}\frac{\rho_\infty}{u_h c_{s,\infty}^3},
155 : * \f}
156 : *
157 : * which gives \cite RezzollaBook
158 : *
159 : * \f{align*}{
160 : * \rho_h\simeq 0.08\frac{\rho_\infty}{c_{s,\infty}^3}.
161 : * \f}
162 : *
163 : * The magnetic field \f$b^2\f$ is the same as the \f$\gamma\ne5/3\f$.
164 : */
165 1 : class BondiMichel : public virtual evolution::initial_data::InitialData,
166 : public AnalyticSolution,
167 : public MarkAsAnalyticSolution {
168 : protected:
169 : template <typename DataType>
170 : struct IntermediateVars;
171 :
172 : public:
173 0 : using equation_of_state_type = EquationsOfState::PolytropicFluid<true>;
174 :
175 : /// The mass of the black hole.
176 1 : struct Mass {
177 0 : using type = double;
178 0 : static constexpr Options::String help = {"Mass of the black hole."};
179 0 : static type lower_bound() { return 0.0; }
180 : };
181 :
182 : /// The radius at which the fluid becomes supersonic.
183 1 : struct SonicRadius {
184 0 : using type = double;
185 0 : static constexpr Options::String help = {
186 : "Schwarzschild radius where fluid becomes supersonic."};
187 0 : static type lower_bound() { return 0.0; }
188 : };
189 :
190 : /// The rest mass density of the fluid at the sonic radius.
191 1 : struct SonicDensity {
192 0 : using type = double;
193 0 : static constexpr Options::String help = {
194 : "The density of the fluid at the sonic radius."};
195 0 : static type lower_bound() { return 0.0; }
196 : };
197 :
198 : /// The polytropic exponent for the polytropic fluid.
199 1 : struct PolytropicExponent {
200 0 : using type = double;
201 0 : static constexpr Options::String help = {
202 : "The polytropic exponent for the polytropic fluid."};
203 0 : static type lower_bound() { return 1.0; }
204 : };
205 :
206 : /// The strength of the radial magnetic field.
207 1 : struct MagFieldStrength {
208 0 : using type = double;
209 0 : static constexpr Options::String help = {
210 : "The strength of the radial magnetic field."};
211 : };
212 :
213 0 : using options = tmpl::list<Mass, SonicRadius, SonicDensity,
214 : PolytropicExponent, MagFieldStrength>;
215 0 : static constexpr Options::String help = {
216 : "Bondi-Michel solution with a radial magnetic field using \n"
217 : "the Schwarzschild coordinate system. Quantities prefixed with \n"
218 : "`sonic` refer to field quantities evaluated at the radius \n"
219 : "where the fluid speed overtakes the sound speed."};
220 :
221 0 : BondiMichel() = default;
222 0 : BondiMichel(const BondiMichel& /*rhs*/) = default;
223 0 : BondiMichel& operator=(const BondiMichel& /*rhs*/) = default;
224 0 : BondiMichel(BondiMichel&& /*rhs*/) = default;
225 0 : BondiMichel& operator=(BondiMichel&& /*rhs*/) = default;
226 0 : ~BondiMichel() override = default;
227 :
228 0 : BondiMichel(double mass, double sonic_radius, double sonic_density,
229 : double polytropic_exponent, double mag_field_strength);
230 :
231 0 : auto get_clone() const
232 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
233 :
234 : /// \cond
235 : explicit BondiMichel(CkMigrateMessage* msg);
236 : using PUP::able::register_constructor;
237 : WRAPPED_PUPable_decl_template(BondiMichel);
238 : /// \endcond
239 :
240 : /// Retrieve a collection of hydro variables at `(x, t)`
241 : template <typename DataType, typename... Tags>
242 1 : tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
243 : const double /*t*/,
244 : tmpl::list<Tags...> /*meta*/) const {
245 : // non-const so we can move out the metric vars. We assume that no variable
246 : // is being retrieved more than once, which would cause problems with
247 : // TaggedTuple anyway.
248 : auto intermediate_vars = IntermediateVars<DataType>{
249 : rest_mass_density_at_infinity_,
250 : mass_accretion_rate_over_four_pi_,
251 : mass_,
252 : polytropic_constant_,
253 : polytropic_exponent_,
254 : bernoulli_constant_squared_minus_one_,
255 : sonic_radius_,
256 : sonic_density_,
257 : x,
258 : tmpl2::flat_any_v<
259 : not tmpl::list_contains_v<hydro::grmhd_tags<DataType>, Tags>...>,
260 : background_spacetime_};
261 : return {get<Tags>(variables(x, tmpl::list<Tags>{}, intermediate_vars))...};
262 : }
263 :
264 : template <typename DataType, typename Tag>
265 0 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
266 : const double /*t*/, // NOLINT
267 : tmpl::list<Tag> /*meta*/) const {
268 : return variables(
269 : x, tmpl::list<Tag>{},
270 : IntermediateVars<DataType>{
271 : rest_mass_density_at_infinity_, mass_accretion_rate_over_four_pi_,
272 : mass_, polytropic_constant_, polytropic_exponent_,
273 : bernoulli_constant_squared_minus_one_, sonic_radius_,
274 : sonic_density_, x,
275 : not tmpl::list_contains_v<hydro::grmhd_tags<DataType>, Tag>,
276 : background_spacetime_});
277 : }
278 :
279 : // NOLINTNEXTLINE(google-runtime-references)
280 0 : void pup(PUP::er& /*p*/) override;
281 0 : const EquationsOfState::PolytropicFluid<true>& equation_of_state() const {
282 : return equation_of_state_;
283 : }
284 :
285 : protected:
286 0 : friend bool operator==(const BondiMichel& lhs, const BondiMichel& rhs);
287 :
288 : template <typename DataType>
289 0 : auto variables(const tnsr::I<DataType, 3>& x,
290 : tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/,
291 : const IntermediateVars<DataType>& vars) const
292 : -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
293 :
294 : template <typename DataType>
295 0 : auto variables(const tnsr::I<DataType, 3>& x,
296 : tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/,
297 : const IntermediateVars<DataType>& vars) const
298 : -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
299 :
300 : template <typename DataType>
301 0 : auto variables(
302 : const tnsr::I<DataType, 3>& x,
303 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/,
304 : const IntermediateVars<DataType>& vars) const
305 : -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
306 :
307 : template <typename DataType>
308 0 : auto variables(const tnsr::I<DataType, 3>& x,
309 : tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/,
310 : const IntermediateVars<DataType>& vars) const
311 : -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
312 :
313 : template <typename DataType>
314 0 : auto variables(const tnsr::I<DataType, 3>& x,
315 : tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/,
316 : const IntermediateVars<DataType>& vars) const
317 : -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>>;
318 :
319 : template <typename DataType>
320 0 : auto variables(const tnsr::I<DataType, 3>& x,
321 : tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/,
322 : const IntermediateVars<DataType>& vars) const
323 : -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
324 :
325 : template <typename DataType>
326 0 : auto variables(const tnsr::I<DataType, 3>& x,
327 : tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/,
328 : const IntermediateVars<DataType>& vars) const
329 : -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
330 :
331 : template <typename DataType>
332 0 : auto variables(
333 : const tnsr::I<DataType, 3>& x,
334 : tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/,
335 : const IntermediateVars<DataType>& vars) const
336 : -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
337 :
338 : template <typename DataType>
339 0 : auto variables(const tnsr::I<DataType, 3>& x,
340 : tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/,
341 : const IntermediateVars<DataType>& vars) const
342 : -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
343 :
344 : template <typename DataType>
345 0 : auto variables(const tnsr::I<DataType, 3>& x,
346 : tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/,
347 : const IntermediateVars<DataType>& vars) const
348 : -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
349 :
350 : template <typename DataType, typename Tag,
351 : Requires<not tmpl::list_contains_v<
352 : tmpl::push_back<hydro::grmhd_tags<DataType>,
353 : hydro::Tags::SpecificEnthalpy<DataType>>,
354 : Tag>> = nullptr>
355 0 : tuples::TaggedTuple<Tag> variables(
356 : const tnsr::I<DataType, 3>& /*x*/, tmpl::list<Tag> /*meta*/,
357 : const IntermediateVars<DataType>& vars) const {
358 : return {std::move(get<Tag>(vars.kerr_schild_soln))};
359 : }
360 :
361 : template <typename DataType>
362 0 : struct IntermediateVars {
363 0 : IntermediateVars(double rest_mass_density_at_infinity,
364 : double in_mass_accretion_rate_over_four_pi, double in_mass,
365 : double in_polytropic_constant,
366 : double in_polytropic_exponent,
367 : double in_bernoulli_constant_squared_minus_one,
368 : double in_sonic_radius, double in_sonic_density,
369 : const tnsr::I<DataType, 3>& x, bool need_spacetime,
370 : const gr::Solutions::KerrSchild& background_spacetime);
371 0 : DataType radius{};
372 0 : DataType rest_mass_density{};
373 0 : DataType electron_fraction{};
374 0 : double mass_accretion_rate_over_four_pi{};
375 0 : double mass{};
376 0 : double polytropic_constant{};
377 0 : double polytropic_exponent{};
378 0 : double bernoulli_constant_squared_minus_one{};
379 0 : double sonic_radius{};
380 0 : double sonic_density{};
381 0 : double bernoulli_root_function(double rest_mass_density_guess,
382 : double current_radius) const;
383 : tuples::tagged_tuple_from_typelist<
384 : typename gr::Solutions::KerrSchild::tags<DataType>>
385 0 : kerr_schild_soln{};
386 : };
387 :
388 0 : double mass_ = std::numeric_limits<double>::signaling_NaN();
389 0 : double sonic_radius_ = std::numeric_limits<double>::signaling_NaN();
390 0 : double sonic_density_ = std::numeric_limits<double>::signaling_NaN();
391 0 : double polytropic_exponent_ = std::numeric_limits<double>::signaling_NaN();
392 0 : double mag_field_strength_ = std::numeric_limits<double>::signaling_NaN();
393 0 : double sonic_fluid_speed_squared_ =
394 : std::numeric_limits<double>::signaling_NaN();
395 0 : double sonic_sound_speed_squared_ =
396 : std::numeric_limits<double>::signaling_NaN();
397 0 : double polytropic_constant_ = std::numeric_limits<double>::signaling_NaN();
398 0 : double mass_accretion_rate_over_four_pi_ =
399 : std::numeric_limits<double>::signaling_NaN();
400 0 : double bernoulli_constant_squared_minus_one_ =
401 : std::numeric_limits<double>::signaling_NaN();
402 0 : double rest_mass_density_at_infinity_ =
403 : std::numeric_limits<double>::signaling_NaN();
404 0 : EquationsOfState::PolytropicFluid<true> equation_of_state_{};
405 0 : gr::Solutions::KerrSchild background_spacetime_{};
406 : };
407 :
408 0 : bool operator!=(const BondiMichel& lhs, const BondiMichel& rhs);
409 :
410 : } // namespace grmhd::Solutions
|
