Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <memory>
7 :
8 : #include "DataStructures/TaggedTuple.hpp"
9 : #include "DataStructures/Tensor/TypeAliases.hpp"
10 : #include "Evolution/Systems/ForceFree/Tags.hpp"
11 : #include "Options/Context.hpp"
12 : #include "Options/String.hpp"
13 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
14 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
15 : #include "PointwiseFunctions/GeneralRelativity/KerrSchildCoords.hpp"
16 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
17 : #include "Utilities/Serialization/CharmPupable.hpp"
18 : #include "Utilities/TMPL.hpp"
19 :
20 : /// \cond
21 : class DataVector;
22 : namespace PUP {
23 : class er;
24 : } // namespace PUP
25 : /// \endcond
26 :
27 : namespace ForceFree::Solutions {
28 :
29 : /*!
30 : * \brief An exact electrovacuum force-free solution of Maxwell's equations in
31 : * the Schwarzschild spacetime by Wald \cite Wald1974.
32 : *
33 : * The solution is given in terms of the electromagnetic 4-potential
34 : *
35 : * \begin{equation}
36 : * A_\mu = \frac{B_0}{2}(\phi_\mu + 2a t_\mu)
37 : * \end{equation}
38 : *
39 : * where $B_0$ is the vector potential amplitude, $\phi^\mu = \partial_\phi$,
40 : * $t^\mu = \partial_t$, and $a$ is the (dimensionless) spin of the black hole.
41 : * The case $a=0$ is force-free outside the horizon.
42 : *
43 : * \note This solution is not force-free inside the horizon; the condition
44 : * $E_iB^i = 0$ is still satisfied, but $B^2 > E^2$ is not.
45 : *
46 : * In the spherical Kerr-Schild coordinates, the only nonzero component of
47 : * vector potential is
48 : *
49 : * \begin{equation}
50 : * A_\phi = \frac{B_0}{2}r^2 \sin^2 \theta.
51 : * \end{equation}
52 : *
53 : * Computing magnetic fields,
54 : * \begin{align}
55 : * \tilde{B}^r & = \partial_\theta A_\phi = B_0 r^2 \sin\theta\cos\theta \\
56 : * \tilde{B}^\theta & = - \partial_r A_\phi = -B_0 r \sin^2 \theta \\
57 : * \tilde{B}^\phi &= 0 ,
58 : * \end{align}
59 : *
60 : * Transformation to the Cartesian coordinates gives
61 : * \begin{equation}
62 : * \tilde{B}^x = 0 , \quad
63 : * \tilde{B}^y = 0 , \quad
64 : * \tilde{B}^z = B_0 .
65 : * \end{equation}
66 : *
67 : * Electric fields are given by
68 : *
69 : * \begin{equation}
70 : * E_i = F_{ia}n^a = \frac{1}{\alpha}(F_{i0} - F_{ij}\beta^j) .
71 : * \end{equation}
72 : *
73 : * We omit the derivation and write out results below:
74 : * \begin{equation}
75 : * \tilde{E}^x = - \frac{2 M B_0 y}{r^2}, \quad
76 : * \tilde{E}^y = \frac{2 M B_0 x}{r^2}, \quad
77 : * \tilde{E}^z = 0
78 : * \end{equation}
79 : *
80 : * Note that $\tilde{B}^i \equiv \sqrt{\gamma}B^i$, $\tilde{E}^i \equiv
81 : * \sqrt{\gamma}E^i$, and $\gamma = 1 + 2M/r$ in the (Cartesian) Kerr-Schild
82 : * coordinates. We use $M=1$ Schwarzschild black hole in the Kerr-Schild
83 : * coordinates (see the documentation of gr::Solutions::KerrSchild).
84 : *
85 : */
86 1 : class ExactWald : public evolution::initial_data::InitialData,
87 : public MarkAsAnalyticSolution {
88 : public:
89 0 : struct MagneticFieldAmplitude {
90 0 : using type = double;
91 0 : static constexpr Options::String help = {
92 : "Amplitude of magnetic field along z axis."};
93 : };
94 :
95 0 : using options = tmpl::list<MagneticFieldAmplitude>;
96 :
97 0 : static constexpr Options::String help{
98 : "Exact vacuum solution of Maxwell's equations in Schwarzschild BH "
99 : "spacetime by Wald (1974)."};
100 :
101 0 : ExactWald() = default;
102 0 : ExactWald(const ExactWald&) = default;
103 0 : ExactWald& operator=(const ExactWald&) = default;
104 0 : ExactWald(ExactWald&&) = default;
105 0 : ExactWald& operator=(ExactWald&&) = default;
106 0 : ~ExactWald() override = default;
107 :
108 0 : explicit ExactWald(double magnetic_field_amplitude);
109 :
110 0 : auto get_clone() const
111 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
112 :
113 : /// \cond
114 : explicit ExactWald(CkMigrateMessage* msg);
115 : using PUP::able::register_constructor;
116 : WRAPPED_PUPable_decl_template(ExactWald);
117 : /// \endcond
118 :
119 : // NOLINTNEXTLINE(google-runtime-references)
120 0 : void pup(PUP::er& p) override;
121 :
122 : /// @{
123 : /// Retrieve the EM variables at (x,t).
124 1 : auto variables(const tnsr::I<DataVector, 3>& x, double t,
125 : tmpl::list<Tags::TildeE> /*meta*/) const
126 : -> tuples::TaggedTuple<Tags::TildeE>;
127 :
128 1 : auto variables(const tnsr::I<DataVector, 3>& x, double t,
129 : tmpl::list<Tags::TildeB> /*meta*/) const
130 : -> tuples::TaggedTuple<Tags::TildeB>;
131 :
132 1 : static auto variables(const tnsr::I<DataVector, 3>& x, double t,
133 : tmpl::list<Tags::TildePsi> /*meta*/)
134 : -> tuples::TaggedTuple<Tags::TildePsi>;
135 :
136 1 : static auto variables(const tnsr::I<DataVector, 3>& x, double t,
137 : tmpl::list<Tags::TildePhi> /*meta*/)
138 : -> tuples::TaggedTuple<Tags::TildePhi>;
139 :
140 1 : static auto variables(const tnsr::I<DataVector, 3>& x, double t,
141 : tmpl::list<Tags::TildeQ> /*meta*/)
142 : -> tuples::TaggedTuple<Tags::TildeQ>;
143 : /// @}
144 :
145 : /// Retrieve a collection of EM variables at `(x, t)`
146 : template <typename... Tags>
147 1 : tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataVector, 3>& x,
148 : const double t,
149 : tmpl::list<Tags...> /*meta*/) const {
150 : static_assert(sizeof...(Tags) > 1,
151 : "The generic template will recurse infinitely if only one "
152 : "tag is being retrieved.");
153 : return {get<Tags>(variables(x, t, tmpl::list<Tags>{}))...};
154 : }
155 :
156 : /// Retrieve the metric variables
157 : template <typename Tag>
158 1 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataVector, 3>& x, double t,
159 : tmpl::list<Tag> /*meta*/) const {
160 : return background_spacetime_.variables(x, t, tmpl::list<Tag>{});
161 : }
162 :
163 : private:
164 0 : gr::Solutions::KerrSchild background_spacetime_{
165 : 1.0, {{0.0, 0.0, 0.0}}, {{0.0, 0.0, 0.0}}};
166 0 : double magnetic_field_amplitude_ =
167 : std::numeric_limits<double>::signaling_NaN();
168 :
169 0 : friend bool operator==(const ExactWald& lhs, const ExactWald& rhs);
170 : };
171 :
172 0 : bool operator!=(const ExactWald& lhs, const ExactWald& rhs);
173 :
174 : } // namespace ForceFree::Solutions
|
