FishboneMoncriefDisk.hpp
1 // 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/DataVector.hpp"
11 #include "Options/Options.hpp"
12 #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
13 #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
14 #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp" // IWYU pragma: keep
15 #include "PointwiseFunctions/Hydro/Tags.hpp"
17 #include "Utilities/Requires.hpp"
18 #include "Utilities/TMPL.hpp"
20 
21 /// \cond
22 namespace PUP {
23 class er; // IWYU pragma: keep
24 } // namespace PUP
25 /// \endcond
26 
27 // IWYU pragma: no_include <pup.h>
28 
29 namespace RelativisticEuler {
30 namespace Solutions {
31 
32 /*!
33  * \brief Fluid disk orbiting a Kerr black hole
34  *
35  * The Fishbone-Moncrief solution to the 3D relativistic Euler system
36  * \cite fishbone1976apj, representing the isentropic flow of a thick fluid disk
37  * orbiting a Kerr black hole. In Boyer-Lindquist coordinates
38  * \f$(t, r, \theta, \phi)\f$, the flow is assumed to be purely toroidal,
39  *
40  * \f{align*}
41  * u^\mu = (u^t, 0, 0, u^\phi),
42  * \f}
43  *
44  * where \f$u^\mu\f$ is the 4-velocity. Then, all the fluid quantities are
45  * assumed to share the same symmetries as those of the background spacetime,
46  * namely they are stationary (independent of \f$t\f$), and axially symmetric
47  * (independent of \f$\phi\f$). Note that \f$r\f$ and \f$\theta\f$ can also be
48  * interpreted as Kerr (a.k.a. "spherical" Kerr-Schild) coordinates
49  * (see gr::KerrSchildCoords), and that the symmetries of the equilibrium ensure
50  * that \f$u^\mu\f$ as defined above is also the 4-velocity in Kerr coordinates.
51  *
52  * Self-gravity is neglected, so that the fluid
53  * variables are determined as functions of the metric. Following the treatment
54  * by Kozlowski et al. \cite Kozlowski1978aa (but using signature +2)
55  * the solution is expressed in terms of the quantities
56  *
57  * \f{align*}
58  * \Omega &= \dfrac{u^\phi}{u^t},\\
59  * W &= W_\text{in} - \int_{p_\text{in}}^p\frac{dp}{e + p},
60  * \f}
61  *
62  * where \f$\Omega\f$ is the angular velocity, \f$p\f$ is the fluid pressure,
63  * \f$e\f$ is the energy density, and \f$W\f$ is an auxiliary quantity
64  * interpreted in the Newtonian limit as the total (gravitational + centrifugal)
65  * potential. \f$W_\text{in}\f$ and \f$p_\text{in}\f$ are the potential and
66  * the pressure at the radius of the inner edge, i.e. the closest edge to the
67  * black hole. Here we assume \f$p_\text{in} = 0.\f$ The solution to the Euler
68  * equation is then
69  *
70  * \f{align*}
71  * (u^t)^2 &= \frac{A}{2\Delta\Sigma}\left(1 +
72  * \sqrt{1 + \frac{4l^2\Delta \Sigma^2}{A^2\sin^2\theta}}\right)\\
73  * \Omega &= \frac{\Sigma}{A (u^t)^2}\frac{l}{\sin^2\theta} + \frac{2Mra}{A}\\
74  * u^\phi &= \Omega u^t\\
75  * W &= l\Omega - \ln u^t,
76  * \f}
77  *
78  * where
79  *
80  * \f{align*}
81  * \Sigma = r^2 + a^2\cos^2\theta\qquad
82  * \Delta = r^2 - 2Mr + a^2\qquad
83  * A = (r^2 + a^2)^2 - \Delta a^2 \sin^2\theta
84  * \f}
85  *
86  * and \f$l = u_\phi u^t\f$ is the so-called angular momentum per unit
87  * intertial mass, which is a parameter defining an
88  * individual disk. In deriving the solution, an integration constant has been
89  * chosen so that \f$ W\longrightarrow 0\f$ as \f$r\longrightarrow \infty\f$,
90  * in accordance with the Newtonian limit. Note that, from its definition,
91  * equipotential contours coincide with isobaric contours. Physically, the
92  * matter can fill each of the closed surfaces \f$W = \text{const}\f$, giving
93  * rise to an orbiting thick disk. For \f$W > 0\f$, all equipotentials are open,
94  * whereas for \f$W < 0\f$, some of them will be closed. Should a disk exist,
95  * the pressure reaches a maximum value on the equator at a coordinate radius
96  * \f$r_\text{max}\f$ that is related to the angular momentum per unit inertial
97  * mass via
98  *
99  * \f{align*}
100  * l = \dfrac{M^{1/2}(r_\text{max}^{3/2} + aM^{1/2})(a^2 - 2aM^{1/2}
101  * r_\text{max}^{1/2} + r_\text{max}^2)}{2aM^{1/2}r_\text{max}^{3/2} +
102  * (r_\text{max} - 3M)r_\text{max}^2}.
103  * \f}
104  *
105  * Once \f$W\f$ is determined, an equation of state is required in order to
106  * obtain the thermodynamic variables. If the flow is isentropic, the specific
107  * enthalpy can readily be obtained from the first and second laws of
108  * thermodynamics: one has
109  *
110  * \f{align*}
111  * \frac{dp}{e + p} = \frac{dh}{h}
112  * \f}
113  *
114  * so that
115  *
116  * \f{align*}
117  * h = h_\text{in}\exp(W_\text{in} - W),
118  * \f}
119  *
120  * and the pressure can be obtained from a thermodynamic relation of the form
121  * \f$h = h(p)\f$. Here we assume a polytropic relation
122  *
123  * \f{align*}
124  * p = K\rho^\gamma.
125  * \f}
126  *
127  * Once all the variables are known in Boyer-Lindquist (or Kerr) coordinates, it
128  * is straightforward to write them in Cartesian Kerr-Schild coordinates. The
129  * coordinate transformation in gr::KerrSchildCoords helps read the Jacobian
130  * matrix, which, applied to the azimuthal flow of the disk, gives
131  *
132  * \f{align*}
133  * u_\text{KS}^\mu = u^t(1, -y\Omega, x\Omega, 0),
134  * \f}
135  *
136  * where \f$u^t\f$ and \f$\Omega\f$ are now understood as functions of the
137  * Kerr-Schild coordinates. Finally, the spatial velocity can be readily
138  * obtained from its definition,
139  *
140  * \f{align*}
141  * \alpha v^i = \frac{u^i}{u^t} + \beta^i,
142  * \f}
143  *
144  * where \f$\alpha\f$ and \f$\beta^i\f$ are the lapse and the shift,
145  * respectively.
146  *
147  * \note Kozlowski et al. \cite Kozlowski1978aa denote
148  * \f$l_* = u_\phi u^t\f$ in order to
149  * distinguish this quantity from their own definition \f$l = - u_\phi/u_t\f$.
150  */
152  template <typename DataType, bool NeedSpacetime>
153  struct IntermediateVariables;
154 
155  public:
157 
158  /// The mass of the black hole, \f$M\f$.
159  struct BhMass {
160  using type = double;
161  static constexpr OptionString help = {"The mass of the black hole."};
162  static type lower_bound() noexcept { return 0.0; }
163  };
164  /// The dimensionless black hole spin, \f$\chi = a/M\f$.
165  struct BhDimlessSpin {
166  using type = double;
167  static constexpr OptionString help = {"The dimensionless black hole spin."};
168  static type lower_bound() { return 0.0; }
169  static type upper_bound() { return 1.0; }
170  };
171  /// The radial coordinate of the inner edge of the disk, in units of \f$M\f$.
173  using type = double;
174  static constexpr OptionString help = {
175  "The radial coordinate of the inner edge of the disk."};
176  };
177  /// The radial coordinate of the maximum pressure, in units of \f$M\f$.
179  using type = double;
180  static constexpr OptionString help = {
181  "The radial coordinate of the maximum pressure."};
182  };
183  /// The polytropic constant of the fluid.
185  using type = double;
186  static constexpr OptionString help = {
187  "The polytropic constant of the fluid."};
188  static type lower_bound() noexcept { return 0.; }
189  };
190  /// The polytropic exponent of the fluid.
192  using type = double;
193  static constexpr OptionString help = {
194  "The polytropic exponent of the fluid."};
195  static type lower_bound() noexcept { return 1.; }
196  };
197 
198  using options =
201  static constexpr OptionString help = {
202  "Fluid disk orbiting a Kerr black hole."};
203 
204  FishboneMoncriefDisk() = default;
205  FishboneMoncriefDisk(const FishboneMoncriefDisk& /*rhs*/) = delete;
206  FishboneMoncriefDisk& operator=(const FishboneMoncriefDisk& /*rhs*/) = delete;
207  FishboneMoncriefDisk(FishboneMoncriefDisk&& /*rhs*/) noexcept = default;
208  FishboneMoncriefDisk& operator=(FishboneMoncriefDisk&& /*rhs*/) noexcept =
209  default;
210  ~FishboneMoncriefDisk() = default;
211 
212  FishboneMoncriefDisk(double bh_mass, double bh_dimless_spin,
213  double inner_edge_radius, double max_pressure_radius,
214  double polytropic_constant,
215  double polytropic_exponent) noexcept;
216 
217  // Eventually, if we implement a gr::Solutions::BoyerLindquist
218  // black hole, the following two functions aren't needed, and the algebra
219  // of the three functions after can be simplified by using the corresponding
220  // lapse, shift, and spatial metric.
221  template <typename DataType>
222  DataType sigma(const DataType& r_sqrd, const DataType& sin_theta_sqrd) const
223  noexcept;
224 
225  template <typename DataType>
226  DataType inv_ucase_a(const DataType& r_sqrd, const DataType& sin_theta_sqrd,
227  const DataType& delta) const noexcept;
228 
229  template <typename DataType>
230  DataType four_velocity_t_sqrd(const DataType& r_sqrd,
231  const DataType& sin_theta_sqrd) const noexcept;
232 
233  template <typename DataType>
234  DataType angular_velocity(const DataType& r_sqrd,
235  const DataType& sin_theta_sqrd) const noexcept;
236 
237  template <typename DataType>
238  DataType potential(const DataType& r_sqrd,
239  const DataType& sin_theta_sqrd) const noexcept;
240 
241  SPECTRE_ALWAYS_INLINE size_t index_helper(tmpl::no_such_type_ /*meta*/) const
242  noexcept {
244  }
245  template <typename T>
246  SPECTRE_ALWAYS_INLINE size_t index_helper(T /*meta*/) const noexcept {
247  return T::value;
248  }
249 
250  // @{
251  /// The fluid variables in Cartesian Kerr-Schild coordinates at `(x, t)`
252  ///
253  /// \note The functions are optimized for retrieving the hydro variables
254  /// before the metric variables.
255  template <typename DataType, typename... Tags>
257  const tnsr::I<DataType, 3>& x,
258  const double t, // NOLINT(readability-avoid-const-params-in-decls)
259  tmpl::list<Tags...> /*meta*/) const noexcept {
260  // Can't store IntermediateVariables as member variable because we need to
261  // be threadsafe.
262  IntermediateVariables<
263  DataType,
265  cpp17::is_same_v<Tags, hydro::Tags::SpatialVelocity<DataType, 3>> or
266  cpp17::is_same_v<Tags, hydro::Tags::LorentzFactor<DataType>> or
267  not tmpl::list_contains_v<hydro::grmhd_tags<DataType>, Tags>)...>>
268  vars(bh_spin_a_, background_spacetime_, x, t,
269  index_helper(
270  tmpl::index_of<tmpl::list<Tags...>,
272  index_helper(
273  tmpl::index_of<tmpl::list<Tags...>,
275  return {std::move(get<Tags>(
276  variables(x, tmpl::list<Tags>{}, vars,
277  tmpl::index_of<tmpl::list<Tags...>, Tags>::value)))...};
278  }
279 
280  template <typename DataType, typename Tag>
281  tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
282  const double t, // NOLINT
283  tmpl::list<Tag> /*meta*/) const noexcept {
284  // Can't store IntermediateVariables as member variable because we need to
285  // be threadsafe.
286  IntermediateVariables<
287  DataType,
288  cpp17::is_same_v<Tag, hydro::Tags::SpatialVelocity<DataType, 3>> or
289  cpp17::is_same_v<Tag, hydro::Tags::LorentzFactor<DataType>> or
290  not tmpl::list_contains_v<hydro::grmhd_tags<DataType>, Tag>>
291  intermediate_vars(bh_spin_a_, background_spacetime_, x, t,
294  return variables(x, tmpl::list<Tag>{}, intermediate_vars, 0);
295  }
296  // @}
297 
298  // clang-tidy: no runtime references
299  void pup(PUP::er& /*p*/) noexcept; // NOLINT
300 
301  const EquationsOfState::PolytropicFluid<true>& equation_of_state() const
302  noexcept {
303  return equation_of_state_;
304  }
305 
306  private:
307  template <typename DataType, bool NeedSpacetime>
308  auto variables(const tnsr::I<DataType, 3>& x,
309  tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/,
310  const IntermediateVariables<DataType, NeedSpacetime>& vars,
311  size_t index) const noexcept
313 
314  template <typename DataType, bool NeedSpacetime>
315  auto variables(const tnsr::I<DataType, 3>& x,
316  tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/,
317  const IntermediateVariables<DataType, NeedSpacetime>& vars,
318  size_t index) const noexcept
320 
321  template <typename DataType, bool NeedSpacetime>
322  auto variables(const tnsr::I<DataType, 3>& x,
323  tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/,
324  const IntermediateVariables<DataType, NeedSpacetime>& vars,
325  size_t index) const noexcept
327 
328  template <typename DataType, bool NeedSpacetime>
329  auto variables(
330  const tnsr::I<DataType, 3>& x,
332  const IntermediateVariables<DataType, NeedSpacetime>& vars,
333  size_t index) const noexcept
335 
336  template <typename DataType>
337  auto variables(const tnsr::I<DataType, 3>& x,
338  tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/,
339  const IntermediateVariables<DataType, true>& vars,
340  size_t index) const noexcept
342 
343  template <typename DataType>
344  auto variables(const tnsr::I<DataType, 3>& x,
345  tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/,
346  const IntermediateVariables<DataType, true>& vars,
347  size_t index) const noexcept
349 
350  template <typename DataType, bool NeedSpacetime>
351  auto variables(
352  const tnsr::I<DataType, 3>& x,
353  tmpl::list<
355  const IntermediateVariables<DataType, NeedSpacetime>& vars,
356  size_t index) const noexcept
359 
360  template <typename DataType, bool NeedSpacetime>
361  auto variables(
362  const tnsr::I<DataType, 3>& x,
364  const IntermediateVariables<DataType, NeedSpacetime>& vars,
365  size_t index) const noexcept
367 
368  // Grab the metric variables
369  template <typename DataType, typename Tag,
371  Tag>> = nullptr>
372  tuples::TaggedTuple<Tag> variables(
373  const tnsr::I<DataType, 3>& /*x*/, tmpl::list<Tag> /*meta*/,
374  IntermediateVariables<DataType, true>& vars, const size_t index) const
375  noexcept {
376  // The only hydro variables that use the GR solution are spatial velocity
377  // and Lorentz factor. If those have already been set (their index in the
378  // typelist is lower than our index) then we can safely std::move the GR
379  // solution out, assuming nobody gets the same variable twice (e.g. trying
380  // to retrieve the lapse twice from the solution in the same call to the
381  // function), but then TaggedTuple would explode horribly, so nothing to
382  // worry about.
383  // Analytic solutions are used in Dirichlet boundary conditions and thus are
384  // called quite frequently, making some optimization worthwhile.
385  if (index > vars.spatial_velocity_index and
386  index > vars.lorentz_factor_index) {
387  return {std::move(get<Tag>(vars.kerr_schild_soln))};
388  }
389  return {get<Tag>(vars.kerr_schild_soln)};
390  }
391 
392  template <typename DataType, bool NeedSpacetime, typename Func>
393  void variables_impl(
394  const IntermediateVariables<DataType, NeedSpacetime>& vars, Func f) const
395  noexcept;
396 
397  // Intermediate variables needed to set several of the Fishbone-Moncrief
398  // solution's variables.
399  template <typename DataType, bool NeedSpacetime>
400  struct IntermediateVariables {
401  IntermediateVariables(double bh_spin_a,
402  const gr::Solutions::KerrSchild& background_spacetime,
403  const tnsr::I<DataType, 3>& x, double t,
404  size_t in_spatial_velocity_index,
405  size_t in_lorentz_factor_index) noexcept;
406 
407  DataType r_squared{};
408  DataType sin_theta_squared{};
409  tuples::tagged_tuple_from_typelist<
410  typename gr::Solutions::KerrSchild::tags<DataType>>
411  kerr_schild_soln{};
412  size_t spatial_velocity_index;
413  size_t lorentz_factor_index;
414  };
415 
416  friend bool operator==(const FishboneMoncriefDisk& lhs,
417  const FishboneMoncriefDisk& rhs) noexcept;
418 
419  double bh_mass_ = std::numeric_limits<double>::signaling_NaN();
420  double bh_spin_a_ = std::numeric_limits<double>::signaling_NaN();
421  double inner_edge_radius_ = std::numeric_limits<double>::signaling_NaN();
422  double max_pressure_radius_ = std::numeric_limits<double>::signaling_NaN();
423  double polytropic_constant_ = std::numeric_limits<double>::signaling_NaN();
424  double polytropic_exponent_ = std::numeric_limits<double>::signaling_NaN();
425  double angular_momentum_ = std::numeric_limits<double>::signaling_NaN();
426  EquationsOfState::PolytropicFluid<true> equation_of_state_{};
427  gr::Solutions::KerrSchild background_spacetime_{};
428 };
429 
430 bool operator!=(const FishboneMoncriefDisk& lhs,
431  const FishboneMoncriefDisk& rhs) noexcept;
432 } // namespace Solutions
433 } // namespace RelativisticEuler
Definition: Strahlkorper.hpp:14
Defines class tuples::TaggedTuple.
The spatial velocity .
Definition: Tags.hpp:144
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)
Items related to evolving the relativistic Euler system.
Definition: Characteristics.hpp:21
The radial coordinate of the maximum pressure, in units of .
Definition: FishboneMoncriefDisk.hpp:178
The dimensionless black hole spin, .
Definition: FishboneMoncriefDisk.hpp:165
Defines classes and functions for making classes creatable from input files.
The radial coordinate of the inner edge of the disk, in units of .
Definition: FishboneMoncriefDisk.hpp:172
Defines the type alias Requires.
The polytropic constant of the fluid.
Definition: FishboneMoncriefDisk.hpp:184
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:27
tuples::TaggedTuple< Tags... > variables(const tnsr::I< DataType, 3 > &x, const double t, tmpl::list< Tags... >) const noexcept
The fluid variables in Cartesian Kerr-Schild coordinates at (x, t)
Definition: FishboneMoncriefDisk.hpp:256
#define SPECTRE_ALWAYS_INLINE
Always inline a function. Only use this if you benchmarked the code.
Definition: ForceInline.hpp:20
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
T max(T... args)
The Lorentz factor .
Definition: Tags.hpp:64
Defines macro to always inline a function.
Definition: DataBoxTag.hpp:29
Defines a list of useful type aliases for tensors.
Fluid disk orbiting a Kerr black hole.
Definition: FishboneMoncriefDisk.hpp:151
Wraps the template metaprogramming library used (brigand)
Kerr black hole in Kerr-Schild coordinates.
Definition: KerrSchild.hpp:209
The polytropic exponent of the fluid.
Definition: FishboneMoncriefDisk.hpp:191
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
tuples::TaggedTuple< Tag > variables(const tnsr::I< DataType, 3 > &x, const double t, tmpl::list< Tag >) const noexcept
The fluid variables in Cartesian Kerr-Schild coordinates at (x, t)
Definition: FishboneMoncriefDisk.hpp:281
The specific enthalpy .
Definition: Tags.hpp:169
The mass of the black hole, .
Definition: FishboneMoncriefDisk.hpp:159