BondiMichel.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <limits>
7 #include <pup.h>
8 
9 #include "DataStructures/DataBox/Prefixes.hpp" // IWYU pragma: keep
10 #include "DataStructures/DataVector.hpp"
12 #include "Options/Options.hpp"
13 #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
14 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" // IWYU pragma: keep
15 #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
16 #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp" // IWYU pragma: keep
17 #include "PointwiseFunctions/Hydro/Tags.hpp" // IWYU pragma: keep
18 #include "Utilities/Requires.hpp"
19 #include "Utilities/TMPL.hpp"
21 
22 namespace grmhd {
23 namespace Solutions {
24 
25 /*!
26  * \brief Bondi-Michel accretion with superposed magnetic field in Schwarzschild
27  * spacetime in Kerr-Schild coordinates.
28  *
29  * An analytic solution to the 3-D GrMhd system. The user specifies the sonic
30  * radius \f$r_c\f$ and sonic rest mass density \f$\rho_c\f$, which are the
31  * radius and rest mass density at the sonic point, the radius at which the
32  * fluid's Eulerian velocity as seen by a distant observer overtakes the local
33  * sound speed \f$c_{s,c}\f$. With a specified polytropic exponent \f$\Gamma\f$,
34  * these quantities can be related to the sound speed at infinity
35  * \f$c_{s,\inf}\f$ using the following relations:
36  *
37  * \f{align*}
38  * c_{s,c}^2 &= \frac{1}{2r_c - 3} \\
39  * c_{s,\inf}^2 &= \Gamma - 1 + (c_{s,c}^2 - \Gamma + 1)\sqrt{1 + 3c_{s,c}^2}
40  * \f}
41  *
42  * In the case of the interstellar medium, the sound
43  * speed is \f$\approx 10^{-4}\f$, which results in a sonic radius of
44  * \f$\approx 10^8 M\f$ for \f$\Gamma \neq 5/3\f$ (Rezzola and Zanotti, 2013).
45  * However, for numerical-testing purposes, it is more common to use a value of
46  * \f$\approx 10 M\f$.
47  *
48  * The density is found via root-finding, through the
49  * Bernoulli equation. As one approaches the sonic radius, a second root makes
50  * an appearance and one must take care to bracket the correct root. This is
51  * done by using the upper bound \f$\frac{\dot{M}}{4\pi}\sqrt{\frac{2}{Mr^3}}\f$
52  *
53  * Additionally specified by the user are the polytropic exponent \f$\Gamma\f$,
54  * and the strength parameter of the magnetic field \f$B\f$.
55  * In Kerr-Schild-Cartesian coordinates \f$(x, y, z)\f$, where
56  * \f$ r = \sqrt{x^2 + y^2 + z^2}\f$, the superposed magnetic field is
57  * \f{align*}
58  * B_x(\vec{x},t) &= \frac{B x}{r^3 \sqrt{1 + 2/r}} \\
59  * B_y(\vec{x},t) &= \frac{B y}{r^3 \sqrt{1 + 2/r}} \\
60  * B_z(\vec{x},t) &= \frac{B z}{r^3 \sqrt{1 + 2/r}}
61  * \f}
62  */
63 class BondiMichel {
64  template <typename DataType>
65  struct IntermediateVars;
66 
67  public:
70 
71  /// The mass of the black hole.
72  struct Mass {
73  using type = double;
74  static constexpr OptionString help = {"Mass of the black hole."};
75  static type lower_bound() noexcept { return 0.0; }
76  };
77 
78  /// The radius at which the fluid becomes supersonic.
79  struct SonicRadius {
80  using type = double;
81  static constexpr OptionString help = {
82  "Schwarzschild radius where fluid becomes supersonic."};
83  static type lower_bound() noexcept { return 0.0; }
84  };
85 
86  /// The rest mass density of the fluid at the sonic radius.
87  struct SonicDensity {
88  using type = double;
89  static constexpr OptionString help = {
90  "The density of the fluid at the sonic radius."};
91  static type lower_bound() noexcept { return 0.0; }
92  };
93 
94  /// The polytropic exponent for the polytropic fluid.
96  using type = double;
97  static constexpr OptionString help = {
98  "The polytropic exponent for the polytropic fluid."};
99  static type lower_bound() noexcept { return 1.0; }
100  };
101 
102  /// The strength of the radial magnetic field.
104  using type = double;
105  static constexpr OptionString help = {
106  "The strength of the radial magnetic field."};
107  };
108 
109  using options = tmpl::list<Mass, SonicRadius, SonicDensity,
111  static constexpr OptionString help = {
112  "Bondi-Michel solution with a radial magnetic field using \n"
113  "the Schwarzschild coordinate system. Quantities prefixed with \n"
114  "`sonic` refer to field quantities evaluated at the radius \n"
115  "where the fluid speed overtakes the sound speed."};
116 
117  BondiMichel() = default;
118  BondiMichel(const BondiMichel& /*rhs*/) = delete;
119  BondiMichel& operator=(const BondiMichel& /*rhs*/) = delete;
120  BondiMichel(BondiMichel&& /*rhs*/) noexcept = default;
121  BondiMichel& operator=(BondiMichel&& /*rhs*/) noexcept = default;
122  ~BondiMichel() = default;
123 
124  BondiMichel(Mass::type mass, SonicRadius::type sonic_radius,
125  SonicDensity::type sonic_density,
126  PolytropicExponent::type polytropic_exponent,
127  MagFieldStrength::type mag_field_strength) noexcept;
128 
129  /// Retrieve a collection of hydro variables at `(x, t)`
130  template <typename DataType, typename... Tags>
131  tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
132  const double /*t*/,
133  tmpl::list<Tags...> /*meta*/) const
134  noexcept {
135  // non-const so we can move out the metric vars. We assume that no variable
136  // is being retrieved more than once, which would cause problems with
137  // TaggedTuple anyway.
138  auto intermediate_vars = IntermediateVars<DataType>{
139  rest_mass_density_at_infinity_,
140  mass_accretion_rate_over_four_pi_,
141  mass_,
142  polytropic_constant_,
143  polytropic_exponent_,
144  bernoulli_constant_squared_minus_one_,
145  sonic_radius_,
146  sonic_density_,
147  x,
149  not tmpl::list_contains_v<hydro::grmhd_tags<DataType>, Tags>...>,
150  background_spacetime_};
151  return {get<Tags>(variables(x, tmpl::list<Tags>{}, intermediate_vars))...};
152  }
153 
154  template <typename DataType, typename Tag>
155  tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
156  const double /*t*/, // NOLINT
157  tmpl::list<Tag> /*meta*/) const noexcept {
158  return variables(
159  x, tmpl::list<Tag>{},
160  IntermediateVars<DataType>{
161  rest_mass_density_at_infinity_, mass_accretion_rate_over_four_pi_,
162  mass_, polytropic_constant_, polytropic_exponent_,
163  bernoulli_constant_squared_minus_one_, sonic_radius_,
164  sonic_density_, x,
165  not tmpl::list_contains_v<hydro::grmhd_tags<DataType>, Tag>,
166  background_spacetime_});
167  }
168 
169  // clang-tidy: no runtime references
170  void pup(PUP::er& /*p*/) noexcept; // NOLINT
171  const EquationsOfState::PolytropicFluid<true>& equation_of_state() const
172  noexcept {
173  return equation_of_state_;
174  }
175 
176  private:
177  friend bool operator==(const BondiMichel& lhs,
178  const BondiMichel& rhs) noexcept;
179 
180  template <typename DataType>
181  auto variables(const tnsr::I<DataType, 3>& x,
182  tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/,
183  const IntermediateVars<DataType>& vars) const noexcept
185 
186  template <typename DataType>
187  auto variables(
188  const tnsr::I<DataType, 3>& x,
190  const IntermediateVars<DataType>& vars) const noexcept
192 
193  template <typename DataType>
194  auto variables(const tnsr::I<DataType, 3>& x,
195  tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/,
196  const IntermediateVars<DataType>& vars) const noexcept
198 
199  template <typename DataType>
200  auto variables(
201  const tnsr::I<DataType, 3>& x,
202  tmpl::list<
204  const IntermediateVars<DataType>& vars) const noexcept
207 
208  template <typename DataType>
209  auto variables(
210  const tnsr::I<DataType, 3>& x,
211  tmpl::list<
213  const IntermediateVars<DataType>& vars) const noexcept
216 
217  template <typename DataType>
218  auto variables(
219  const tnsr::I<DataType, 3>& x,
221  const IntermediateVars<DataType>& vars) const noexcept
223 
224  template <typename DataType>
225  auto variables(const tnsr::I<DataType, 3>& x,
226  tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/,
227  const IntermediateVars<DataType>& vars) const noexcept
229 
230  template <typename DataType>
231  auto variables(const tnsr::I<DataType, 3>& x,
232  tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/,
233  const IntermediateVars<DataType>& vars) const noexcept
235 
236  template <typename DataType, typename Tag,
238  Tag>> = nullptr>
239  tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& /*x*/,
240  tmpl::list<Tag> /*meta*/,
241  IntermediateVars<DataType>& vars) const
242  noexcept {
243  return {std::move(get<Tag>(vars.kerr_schild_soln))};
244  }
245 
246  template <typename DataType>
247  struct IntermediateVars {
248  IntermediateVars(
249  double rest_mass_density_at_infinity,
250  double in_mass_accretion_rate_over_four_pi, double in_mass,
251  double in_polytropic_constant, double in_polytropic_exponent,
252  double in_bernoulli_constant_squared_minus_one, double in_sonic_radius,
253  double in_sonic_density, const tnsr::I<DataType, 3>& x,
254  bool need_spacetime,
255  const gr::Solutions::KerrSchild& background_spacetime) noexcept;
256  DataType radius{};
257  DataType rest_mass_density{};
258  double mass_accretion_rate_over_four_pi{};
259  double mass{};
260  double polytropic_constant{};
261  double polytropic_exponent{};
262  double bernoulli_constant_squared_minus_one{};
263  double sonic_radius{};
264  double sonic_density{};
265  double bernoulli_root_function(double rest_mass_density_guess,
266  double current_radius) const noexcept;
267  tuples::tagged_tuple_from_typelist<
268  typename gr::Solutions::KerrSchild::tags<DataType>>
269  kerr_schild_soln{};
270  };
271 
272  Mass::type mass_ = std::numeric_limits<double>::signaling_NaN();
273  SonicRadius::type sonic_radius_ =
275  SonicDensity::type sonic_density_ =
277  PolytropicExponent::type polytropic_exponent_ =
279  MagFieldStrength::type mag_field_strength_ =
281  double sonic_fluid_speed_squared_ =
283  double sonic_sound_speed_squared_ =
285  double polytropic_constant_ = std::numeric_limits<double>::signaling_NaN();
286  double mass_accretion_rate_over_four_pi_ =
288  double bernoulli_constant_squared_minus_one_ =
290  double rest_mass_density_at_infinity_ =
292  EquationsOfState::PolytropicFluid<true> equation_of_state_{};
293  gr::Solutions::KerrSchild background_spacetime_{};
294 };
295 
296 bool operator!=(const BondiMichel& lhs, const BondiMichel& rhs) noexcept;
297 
298 } // namespace Solutions
299 } // namespace grmhd
The radius at which the fluid becomes supersonic.
Definition: BondiMichel.hpp:79
Defines class tuples::TaggedTuple.
The strength of the radial magnetic field.
Definition: BondiMichel.hpp:103
The spatial velocity .
Definition: Tags.hpp:144
The polytropic exponent for the polytropic fluid.
Definition: BondiMichel.hpp:95
Bondi-Michel accretion with superposed magnetic field in Schwarzschild spacetime in Kerr-Schild coord...
Definition: BondiMichel.hpp:63
The fluid pressure .
Definition: Tags.hpp:123
constexpr bool flat_any_v
A non-short-circuiting logical OR between bools &#39;B"".
Definition: TMPL.hpp:528
The specific internal energy .
Definition: Tags.hpp:176
T signaling_NaN(T... args)
Defines classes and functions for making classes creatable from input files.
Define prefixes for DataBox tags.
The rest mass density of the fluid at the sonic radius.
Definition: BondiMichel.hpp:87
Defines the type alias Requires.
The magnetic field measured by an Eulerian observer, where is the normal to the spatial hypersurfac...
Definition: Tags.hpp:80
The divergence-cleaning field .
Definition: Tags.hpp:47
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:26
tuples::TaggedTuple< Tags... > variables(const tnsr::I< DataType, 3 > &x, const double, tmpl::list< Tags... >) const noexcept
Retrieve a collection of hydro variables at (x, t)
Definition: BondiMichel.hpp:131
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
The Lorentz factor .
Definition: Tags.hpp:64
Definition: DataBoxTag.hpp:29
Defines a list of useful type aliases for tensors.
Wraps the template metaprogramming library used (brigand)
Kerr black hole in Kerr-Schild coordinates.
Definition: KerrSchild.hpp:209
The rest-mass density .
Definition: Tags.hpp:130
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t ...
Definition: Requires.hpp:67
The mass of the black hole.
Definition: BondiMichel.hpp:72
The specific enthalpy .
Definition: Tags.hpp:169
Items related to general relativistic magnetohydrodynamics (GRMHD)
Definition: Characteristics.hpp:34