Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <limits>
7 :
8 : #include "DataStructures/Tensor/TypeAliases.hpp"
9 : #include "Options/String.hpp"
10 : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
11 : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
12 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
13 : #include "PointwiseFunctions/GeneralRelativity/KerrSchildCoords.hpp"
14 : #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp" // IWYU pragma: keep
15 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
16 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
17 : #include "Utilities/Requires.hpp"
18 : #include "Utilities/Serialization/CharmPupable.hpp"
19 : #include "Utilities/TMPL.hpp"
20 : #include "Utilities/TaggedTuple.hpp"
21 :
22 : // IWYU pragma: no_include <pup.h>
23 :
24 : /// \cond
25 : namespace PUP {
26 : class er; // IWYU pragma: keep
27 : } // namespace PUP
28 : /// \endcond
29 :
30 : namespace grmhd::AnalyticData {
31 :
32 : /*!
33 : * \brief Analytic initial data for axially symmetric Bondi-Hoyle accretion.
34 : *
35 : * In the context of studying Bondi-Hoyle accretion, i.e. non-spherical
36 : * accretion on to a Kerr black hole moving relative to a gas cloud, this class
37 : * implements the method proposed by \cite Font1998 to initialize the GRMHD
38 : * variables. The fluid quantities are initialized with their (constant) values
39 : * far from the black hole, e.g. \f$\rho = \rho_\infty\f$. Here we assume a
40 : * polytropic equation of state, so only the rest mass density, as well as the
41 : * polytropic constant and the polytropic exponent, are provided as inputs.
42 : * The spatial velocity is initialized using a field that ensures that the
43 : * injected gas reproduces a continuous parallel wind at large distances.
44 : * The direction of this flow is chosen to be along the black hole spin.
45 : * In Kerr (or "spherical Kerr-Schild", see gr::KerrSchildCoords) coordinates,
46 : *
47 : * \f{align*}
48 : * v^r &= \frac{1}{\sqrt{\gamma_{rr}}}v_\infty \cos\theta\\
49 : * v^\theta &= -\frac{1}{\sqrt{\gamma_{\theta\theta}}}v_\infty \sin\theta\\
50 : * v^\phi &= 0.
51 : * \f}
52 : *
53 : * where \f$\gamma_{ij} = g_{ij}\f$ is the spatial metric, and \f$v_\infty\f$
54 : * is the flow speed far from the black hole. Note that
55 : * \f$v_\infty^2 = v_i v^i\f$. Finally, following the work by \cite Penner2011,
56 : * the magnetic field is initialized using Wald's solution to Maxwell's
57 : * equations in Kerr black hole spacetime. In Kerr ("spherical
58 : * Kerr-Schild") coordinates, the spatial components of the Faraday tensor read
59 : *
60 : * \f{align*}
61 : * F_{r\theta} &= a B_0 \left[1 + \frac{2Mr}{\Sigma^2}(r^2 - a^2)\right]
62 : * \sin\theta\cos\theta\\
63 : * F_{\theta\phi} &= B_0\left[\Delta +
64 : * \frac{2Mr}{\Sigma^2}(r^4 - a^4)\right]\sin\theta\cos\theta\\
65 : * F_{\phi r} &= - B_0\left[r + \frac{M a^2}{\Sigma^2}
66 : * (r^2 - a^2\cos^2\theta)(1 + \cos^2\theta)\right]\sin^2\theta.
67 : * \f}
68 : *
69 : * where \f$\Sigma = r^2 + a^2\cos^2\theta\f$ and \f$\Delta = r^2 - 2Mr +
70 : * a^2\f$. The associated Eulerian magnetic field is
71 : *
72 : * \f{align*}
73 : * B^r = \frac{F_{\theta\phi}}{\sqrt\gamma},\quad
74 : * B^\theta = \frac{F_{\phi r}}{\sqrt\gamma},\quad
75 : * B^\phi = \frac{F_{r\theta}}{\sqrt\gamma}.
76 : * \f}
77 : *
78 : * where \f$\gamma = \text{det}(\gamma_{ij})\f$. Wald's solution reproduces a
79 : * uniform magnetic field far from the black hole.
80 : */
81 1 : class BondiHoyleAccretion : public virtual evolution::initial_data::InitialData,
82 : public MarkAsAnalyticData,
83 : public AnalyticDataBase {
84 : public:
85 0 : using equation_of_state_type = EquationsOfState::PolytropicFluid<true>;
86 :
87 : /// The mass of the black hole, \f$M\f$.
88 1 : struct BhMass {
89 0 : using type = double;
90 0 : static constexpr Options::String help = {"The mass of the black hole."};
91 0 : static type lower_bound() { return 0.0; }
92 : };
93 : /// The dimensionless black hole spin, \f$a_* = a/M\f$.
94 1 : struct BhDimlessSpin {
95 0 : using type = double;
96 0 : static constexpr Options::String help = {
97 : "The dimensionless black hole spin."};
98 0 : static type lower_bound() { return -1.0; }
99 0 : static type upper_bound() { return 1.0; }
100 : };
101 : /// The rest mass density of the fluid far from the black hole.
102 1 : struct RestMassDensity {
103 0 : using type = double;
104 0 : static constexpr Options::String help = {
105 : "The asymptotic rest mass density."};
106 0 : static type lower_bound() { return 0.0; }
107 : };
108 : /// The magnitude of the spatial velocity far from the black hole.
109 1 : struct FlowSpeed {
110 0 : using type = double;
111 0 : static constexpr Options::String help = {
112 : "The magnitude of the asymptotic flow velocity."};
113 : };
114 : /// The strength of the magnetic field.
115 1 : struct MagFieldStrength {
116 0 : using type = double;
117 0 : static constexpr Options::String help = {
118 : "The strength of the magnetic field."};
119 : };
120 : /// The polytropic constant of the fluid.
121 1 : struct PolytropicConstant {
122 0 : using type = double;
123 0 : static constexpr Options::String help = {
124 : "The polytropic constant of the fluid."};
125 0 : static type lower_bound() { return 0.0; }
126 : };
127 : /// The polytropic exponent of the fluid.
128 1 : struct PolytropicExponent {
129 0 : using type = double;
130 0 : static constexpr Options::String help = {
131 : "The polytropic exponent of the fluid."};
132 0 : static type lower_bound() { return 1.0; }
133 : };
134 :
135 0 : using options =
136 : tmpl::list<BhMass, BhDimlessSpin, RestMassDensity, FlowSpeed,
137 : MagFieldStrength, PolytropicConstant, PolytropicExponent>;
138 :
139 0 : static constexpr Options::String help = {
140 : "Axially symmetric accretion on to a Kerr black hole."};
141 :
142 0 : BondiHoyleAccretion() = default;
143 0 : BondiHoyleAccretion(const BondiHoyleAccretion& /*rhs*/) = default;
144 0 : BondiHoyleAccretion& operator=(const BondiHoyleAccretion& /*rhs*/) = default;
145 0 : BondiHoyleAccretion(BondiHoyleAccretion&& /*rhs*/) = default;
146 0 : BondiHoyleAccretion& operator=(BondiHoyleAccretion&& /*rhs*/) = default;
147 0 : ~BondiHoyleAccretion() override = default;
148 :
149 0 : BondiHoyleAccretion(double bh_mass, double bh_dimless_spin,
150 : double rest_mass_density, double flow_speed,
151 : double magnetic_field_strength,
152 : double polytropic_constant, double polytropic_exponent);
153 :
154 0 : auto get_clone() const
155 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
156 :
157 : /// \cond
158 : explicit BondiHoyleAccretion(CkMigrateMessage* msg);
159 : using PUP::able::register_constructor;
160 : WRAPPED_PUPable_decl_template(BondiHoyleAccretion);
161 : /// \endcond
162 :
163 : /// @{
164 : /// Retrieve hydro variable at `x`
165 : template <typename DataType>
166 1 : auto variables(const tnsr::I<DataType, 3>& x,
167 : tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/)
168 : const -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
169 :
170 : template <typename DataType>
171 1 : auto variables(const tnsr::I<DataType, 3>& x,
172 : tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/)
173 : const -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
174 :
175 : template <typename DataType>
176 1 : auto variables(
177 : const tnsr::I<DataType, 3>& x,
178 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/) const
179 : -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
180 :
181 : template <typename DataType>
182 1 : auto variables(const tnsr::I<DataType, 3>& x,
183 : tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/) const
184 : -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>>;
185 :
186 : template <typename DataType>
187 1 : auto variables(const tnsr::I<DataType, 3>& x,
188 : tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/) const
189 : -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
190 :
191 : template <typename DataType>
192 1 : auto variables(const tnsr::I<DataType, 3>& x,
193 : tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/)
194 : const -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
195 :
196 : template <typename DataType>
197 1 : auto variables(const tnsr::I<DataType, 3>& x,
198 : tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/)
199 : const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
200 :
201 : template <typename DataType>
202 1 : auto variables(
203 : const tnsr::I<DataType, 3>& x,
204 : tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/) const
205 : -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
206 :
207 : template <typename DataType>
208 1 : auto variables(const tnsr::I<DataType, 3>& x,
209 : tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/)
210 : const -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
211 :
212 : template <typename DataType>
213 1 : auto variables(const tnsr::I<DataType, 3>& x,
214 : tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/)
215 : const -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
216 : /// @}
217 :
218 : /// Retrieve a collection of hydro variables at `x`
219 : template <typename DataType, typename... Tags>
220 1 : tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
221 : tmpl::list<Tags...> /*meta*/) const {
222 : static_assert(sizeof...(Tags) > 1,
223 : "The generic template will recurse infinitely if only one "
224 : "tag is being retrieved.");
225 : return {tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
226 : }
227 :
228 : /// Retrieve the metric variables at `x`
229 : template <typename DataType, typename Tag,
230 : Requires<not tmpl::list_contains_v<
231 : tmpl::push_back<hydro::grmhd_tags<DataType>,
232 : hydro::Tags::SpecificEnthalpy<DataType>>,
233 : Tag>> = nullptr>
234 1 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
235 : tmpl::list<Tag> /*meta*/) const {
236 : constexpr double dummy_time = 0.0;
237 : return {std::move(get<Tag>(background_spacetime_.variables(
238 : x, dummy_time, gr::Solutions::KerrSchild::tags<DataType>{})))};
239 : }
240 :
241 : // NOLINTNEXTLINE(google-runtime-references)
242 0 : void pup(PUP::er& /*p*/) override;
243 :
244 0 : const EquationsOfState::PolytropicFluid<true>& equation_of_state() const {
245 : return equation_of_state_;
246 : }
247 :
248 : private:
249 0 : friend bool operator==(const BondiHoyleAccretion& lhs,
250 : const BondiHoyleAccretion& rhs);
251 :
252 : // compute the spatial velocity in spherical Kerr-Schild coordinates
253 : template <typename DataType>
254 0 : tnsr::I<DataType, 3, Frame::NoFrame> spatial_velocity(
255 : const DataType& r_squared, const DataType& cos_theta,
256 : const DataType& sin_theta) const;
257 : // compute the magnetic field in spherical Kerr-Schild coordinates
258 : template <typename DataType>
259 0 : tnsr::I<DataType, 3, Frame::NoFrame> magnetic_field(
260 : const DataType& r_squared, const DataType& cos_theta,
261 : const DataType& sin_theta) const;
262 :
263 0 : double bh_mass_ = std::numeric_limits<double>::signaling_NaN();
264 0 : double bh_spin_a_ = std::numeric_limits<double>::signaling_NaN();
265 0 : double rest_mass_density_ = std::numeric_limits<double>::signaling_NaN();
266 0 : double flow_speed_ = std::numeric_limits<double>::signaling_NaN();
267 0 : double magnetic_field_strength_ =
268 : std::numeric_limits<double>::signaling_NaN();
269 0 : double polytropic_constant_ = std::numeric_limits<double>::signaling_NaN();
270 0 : double polytropic_exponent_ = std::numeric_limits<double>::signaling_NaN();
271 0 : EquationsOfState::PolytropicFluid<true> equation_of_state_{};
272 0 : gr::Solutions::KerrSchild background_spacetime_{};
273 0 : gr::KerrSchildCoords kerr_schild_coords_{};
274 : };
275 :
276 0 : bool operator!=(const BondiHoyleAccretion& lhs, const BondiHoyleAccretion& rhs);
277 :
278 : } // namespace grmhd::AnalyticData
|
