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