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 : #include <memory>
8 : #include <pup.h>
9 :
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "Evolution/Systems/ForceFree/Tags.hpp"
12 : #include "Options/Options.hpp"
13 : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
14 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp"
15 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
16 : #include "Utilities/TMPL.hpp"
17 :
18 : /// \cond
19 : class DataVector;
20 : namespace PUP {
21 : class er;
22 : } // namespace PUP
23 : /// \endcond
24 :
25 : namespace ForceFree::AnalyticData {
26 : /*!
27 : * \brief The magnetospheric Wald problem proposed in \cite Komissarov2004
28 : *
29 : * This is an initial value problem that evolves the magnetosphere of a rotating
30 : * black hole. The initial condition is given as same as the exact Wald solution
31 : * \cite Wald1974 (see also documentation of ForceFree::Solutions::ExactWald)
32 : *
33 : * \begin{equation}
34 : * A_\mu = \frac{B_0}{2}(\phi_\mu + 2a t_\mu) ,
35 : * \end{equation}
36 : *
37 : * but electric field is set to zero at $t=0$.
38 : *
39 : * In the cartesian projection of the spherical Kerr-Schild coordinates
40 : * (which we use in the code for representing tensors), initial magnetic fields
41 : * is given as
42 : *
43 : * \begin{align}
44 : * \tilde{B}^{x} &= a B_0 z \left[ (ax-ry) \left\{
45 : * \frac{1}{r^4} + \frac{2M r (r^2-a^2)}{(r^4+a^2z^2)^2} \right\}
46 : * + a M r x \left\{ \frac{r^2-z^2}{r^4(r^4+a^2z^2)}
47 : * - \frac{4(r^2+z^2)}{(r^4+a^2z^2)^2} \right\} \right] \\
48 : * \tilde{B}^{y} &= a B_0 z \left[ (rx+ay) \left\{
49 : * \frac{1}{r^4} + \frac{2M r (r^2-a^2)}{(r^4+a^2z^2)^2} \right\}
50 : * + a M r y \left\{ \frac{r^2-z^2}{r^4(r^4+a^2z^2)}
51 : * - \frac{4(r^2+z^2)}{(r^4+a^2z^2)^2} \right\} \right] \\
52 : * \tilde{B}^{z} &= B_0 \left[
53 : * 1 + \frac{a^2z^2}{r^4} + \frac{M a^2}{r^3}\left\{
54 : * 1 - \frac{z^2(a^2+z^2)(5r^4+a^2z^2)}{(r^4+a^2z^2)^2} \right\} \right] .
55 : * \end{align}
56 : *
57 : * where $M$ and $a$ are mass and (dimensionless) spin of the Kerr black hole,
58 : * $B_0$ is the amplitude of magnetic field, and $r$ is the radial coordinate
59 : * defined in the spherical Kerr-Schild coordinate system (see the documentation
60 : * of gr::Solutions::SphericalKerrSchild). All other variables are set to zero
61 : * at $t=0$.
62 : *
63 : * There is no known exact solution to this problem, but numerical simulations
64 : * \cite Komissarov2004 \cite Paschalidis2013 \cite Etienne2017 report that the
65 : * system converges to a steady state with an equatorial current sheet inside
66 : * the ergosphere.
67 : *
68 : * \note We set $M=1$ and $B_0=1$ in the initial data to fix scales.
69 : *
70 : */
71 1 : class MagnetosphericWald : public evolution::initial_data::InitialData,
72 : public MarkAsAnalyticData {
73 : public:
74 0 : struct Spin {
75 0 : using type = double;
76 0 : static constexpr Options::String help = {
77 : "The dimensionless spin of the Kerr BH"};
78 0 : static type upper_bound() { return 1.0; }
79 0 : static type lower_bound() { return -1.0; }
80 : };
81 :
82 0 : using options = tmpl::list<Spin>;
83 0 : static constexpr Options::String help{
84 : "Magnetospheric Wald initial value problem"};
85 :
86 0 : MagnetosphericWald() = default;
87 0 : MagnetosphericWald(const MagnetosphericWald&) = default;
88 0 : MagnetosphericWald& operator=(const MagnetosphericWald&) = default;
89 0 : MagnetosphericWald(MagnetosphericWald&&) = default;
90 0 : MagnetosphericWald& operator=(MagnetosphericWald&&) = default;
91 0 : ~MagnetosphericWald() override = default;
92 :
93 0 : explicit MagnetosphericWald(double spin,
94 : const Options::Context& context = {});
95 :
96 0 : auto get_clone() const
97 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
98 :
99 : /// \cond
100 : explicit MagnetosphericWald(CkMigrateMessage* msg);
101 : using PUP::able::register_constructor;
102 : WRAPPED_PUPable_decl_template(MagnetosphericWald);
103 : /// \endcond
104 :
105 : // NOLINTNEXTLINE(google-runtime-references)
106 0 : void pup(PUP::er& p) override;
107 :
108 : /// @{
109 : /// Retrieve the EM variables at (x,t).
110 1 : static auto variables(const tnsr::I<DataVector, 3>& x,
111 : tmpl::list<Tags::TildeE> /*meta*/)
112 : -> tuples::TaggedTuple<Tags::TildeE>;
113 :
114 1 : auto variables(const tnsr::I<DataVector, 3>& x,
115 : tmpl::list<Tags::TildeB> /*meta*/) const
116 : -> tuples::TaggedTuple<Tags::TildeB>;
117 :
118 1 : static auto variables(const tnsr::I<DataVector, 3>& x,
119 : tmpl::list<Tags::TildePsi> /*meta*/)
120 : -> tuples::TaggedTuple<Tags::TildePsi>;
121 :
122 1 : static auto variables(const tnsr::I<DataVector, 3>& x,
123 : tmpl::list<Tags::TildePhi> /*meta*/)
124 : -> tuples::TaggedTuple<Tags::TildePhi>;
125 :
126 1 : static auto variables(const tnsr::I<DataVector, 3>& x,
127 : tmpl::list<Tags::TildeQ> /*meta*/)
128 : -> tuples::TaggedTuple<Tags::TildeQ>;
129 : /// @}
130 :
131 : /// Retrieve a collection of EM variables at position x
132 : template <typename... Tags>
133 1 : tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataVector, 3>& x,
134 : tmpl::list<Tags...> /*meta*/) const {
135 : static_assert(sizeof...(Tags) > 1,
136 : "The generic template will recurse infinitely if only one "
137 : "tag is being retrieved.");
138 : return {get<Tags>(variables(x, tmpl::list<Tags>{}))...};
139 : }
140 :
141 : /// Retrieve the metric variables
142 : template <typename Tag>
143 1 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataVector, 3>& x,
144 : tmpl::list<Tag> /*meta*/) const {
145 : constexpr double dummy_time = 0.0;
146 : return background_spacetime_.variables(x, dummy_time, tmpl::list<Tag>{});
147 : }
148 :
149 : private:
150 0 : double spin_ = std::numeric_limits<double>::signaling_NaN();
151 0 : gr::Solutions::SphericalKerrSchild background_spacetime_{};
152 :
153 0 : friend bool operator==(const MagnetosphericWald& lhs,
154 : const MagnetosphericWald& rhs);
155 : };
156 :
157 0 : bool operator!=(const MagnetosphericWald& lhs, const MagnetosphericWald& rhs);
158 :
159 : } // namespace ForceFree::AnalyticData
|
