Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <cstdio>
9 : #include <limits>
10 : #include <optional>
11 : #include <pup.h>
12 : #include <vector>
13 :
14 : #include "DataStructures/DataBox/Prefixes.hpp"
15 : #include "DataStructures/TaggedTuple.hpp"
16 : #include "DataStructures/Tensor/Tensor.hpp"
17 : #include "Elliptic/Systems/SelfForce/GeneralRelativity/Tags.hpp"
18 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
19 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
20 : #include "Options/Auto.hpp"
21 : #include "Options/String.hpp"
22 : #include "PointwiseFunctions/InitialDataUtilities/Background.hpp"
23 : #include "PointwiseFunctions/InitialDataUtilities/InitialGuess.hpp"
24 : #include "Utilities/Gsl.hpp"
25 : #include "Utilities/TMPL.hpp"
26 :
27 0 : namespace GrSelfForce::AnalyticData {
28 :
29 : /*!
30 : * \brief Gravitational self-force of a point mass on a circular equatorial
31 : * orbit in Kerr.
32 : *
33 : * This class defines the gravitational self-force equations for a circular
34 : * orbit by setting the coefficients $\alpha$, $\beta$, and $\gamma$ (see
35 : * `GrSelfForce::FirstOrderSystem`). It also sets the effective source
36 : * $S_m^\mathrm{eff}$ and the singular field $\Psi_m^P$ in the regularized
37 : * region. The coefficients are computed using Mathematica-generated functions
38 : * (see CircularOrbitCoeffs.hpp) and the effective source is computed using the
39 : * GravitationalEffectiveSource code by Wardell et. al.
40 : * (https://github.com/barrywardell/GravitationalEffectiveSource) and then
41 : * transformed to our form of the equations with more Mathematica-generated
42 : * functions (see CircularOrbitConvertEffsource.hpp). The derivation of these
43 : * equations will be presented in a future publication. A very strong test of
44 : * the validity of these equations is evaluating them on the singular field and
45 : * the corresponding effective source provided by the external
46 : * GravitationalEffectiveSource code (see Test_CircularOrbit.cpp).
47 : */
48 1 : class CircularOrbit : public elliptic::analytic_data::Background,
49 : public elliptic::analytic_data::InitialGuess {
50 : public:
51 0 : struct BlackHoleMass {
52 0 : static constexpr Options::String help =
53 : "Kerr mass parameter 'M' of the black hole";
54 0 : using type = double;
55 : };
56 0 : struct BlackHoleSpin {
57 0 : static constexpr Options::String help =
58 : "Kerr dimensionless spin parameter 'chi' of the black hole";
59 0 : using type = double;
60 : };
61 0 : struct OrbitalRadius {
62 0 : static constexpr Options::String help =
63 : "Radius 'r_0' of the circular orbit";
64 0 : using type = double;
65 : };
66 0 : struct MModeNumber {
67 0 : static constexpr Options::String help =
68 : "Mode number 'm' of the m-mode decomposition";
69 0 : using type = int;
70 : };
71 0 : struct HyperboloidalSlicingTransitions {
72 0 : static constexpr Options::String help =
73 : "Enable hyperboloidal slicing by specifying the transition points for "
74 : "the boost function. The boost function transitions from -1 to zero "
75 : "between the first two points and from zero to 1 between the last "
76 : "two points. The two points can be the same, in which case the "
77 : "transition "
78 : "is discontinuous (vtu slicing) and the jump must be handled by the DG "
79 : "scheme (see option 'NullSlicingBlocks'). "
80 : "The effective source can only be evaluated where the "
81 : "boost function is zero, so the regularized region must be between "
82 : "the second and third points.";
83 0 : using type = Options::Auto<std::array<double, 4>, Options::AutoLabel::None>;
84 : };
85 0 : struct PenetratingHorizon {
86 0 : static constexpr Options::String help =
87 : "If 'False', use tortoise radial coordinate 'r_star' where the Kerr "
88 : "horizon is at negative infinity, and the angular coordinate is "
89 : "'theta'. If 'True', use Boyer-Lindquist radial coordinate 'r' "
90 : "where the Kerr horizon is at r_+, and the angular coordinate is "
91 : "'cos(theta)'.";
92 0 : using type = bool;
93 : };
94 0 : using options =
95 : tmpl::list<BlackHoleMass, BlackHoleSpin, OrbitalRadius, MModeNumber,
96 : HyperboloidalSlicingTransitions, PenetratingHorizon>;
97 0 : static constexpr Options::String help =
98 : "Quasicircular orbit of a point mass in Kerr spacetime";
99 :
100 0 : CircularOrbit() = default;
101 0 : CircularOrbit(const CircularOrbit&) = default;
102 0 : CircularOrbit& operator=(const CircularOrbit&) = default;
103 0 : CircularOrbit(CircularOrbit&&) = default;
104 0 : CircularOrbit& operator=(CircularOrbit&&) = default;
105 0 : ~CircularOrbit() override = default;
106 :
107 0 : CircularOrbit(double black_hole_mass, double black_hole_spin,
108 : double orbital_radius, int m_mode_number,
109 : const std::optional<std::array<double, 4>>&
110 : hyperboloidal_slicing_transitions,
111 : bool penetrating_horizon);
112 :
113 0 : explicit CircularOrbit(CkMigrateMessage* m);
114 : using PUP::able::register_constructor;
115 0 : WRAPPED_PUPable_decl_template(CircularOrbit);
116 :
117 0 : tnsr::I<double, 2> puncture_position() const;
118 0 : double black_hole_mass() const { return black_hole_mass_; }
119 0 : double black_hole_spin() const { return black_hole_spin_; }
120 0 : double orbital_radius() const { return orbital_radius_; }
121 0 : double omega() const;
122 0 : int m_mode_number() const { return m_mode_number_; }
123 : const std::optional<std::array<double, 4>>&
124 0 : hyperboloidal_slicing_transitions() const {
125 : return hyperboloidal_slicing_transitions_;
126 : }
127 0 : bool penetrating_horizon() const { return penetrating_horizon_; }
128 :
129 0 : using background_tags =
130 : tmpl::list<Tags::Alpha, Tags::Beta, Tags::GammaRstar, Tags::GammaTheta>;
131 0 : using source_tags = tmpl::list<
132 : ::Tags::FixedSource<Tags::MMode>, Tags::SingularField,
133 : ::Tags::deriv<Tags::SingularField, tmpl::size_t<2>, Frame::Inertial>,
134 : Tags::BoyerLindquistRadius>;
135 :
136 : // Background
137 0 : tuples::tagged_tuple_from_typelist<background_tags> variables(
138 : const tnsr::I<DataVector, 2>& x, background_tags /*meta*/) const;
139 :
140 : // Initial guess
141 0 : static tuples::TaggedTuple<Tags::MMode> variables(
142 : const tnsr::I<DataVector, 2>& x, tmpl::list<Tags::MMode> /*meta*/);
143 :
144 : // Fixed sources
145 0 : tuples::tagged_tuple_from_typelist<source_tags> variables(
146 : const tnsr::I<DataVector, 2>& x, source_tags /*meta*/) const;
147 :
148 : template <typename... RequestedTags>
149 0 : tuples::TaggedTuple<RequestedTags...> variables(
150 : const tnsr::I<DataVector, 2>& x, const Mesh<2>& /*mesh*/,
151 : const InverseJacobian<DataVector, 2, Frame::ElementLogical,
152 : Frame::Inertial>& /*inv_jacobian*/,
153 : tmpl::list<RequestedTags...> /*meta*/) const {
154 : return variables(x, tmpl::list<RequestedTags...>{});
155 : }
156 :
157 : // NOLINTNEXTLINE
158 0 : void pup(PUP::er& p) override;
159 :
160 : private:
161 0 : friend bool operator==(const CircularOrbit& lhs, const CircularOrbit& rhs);
162 :
163 0 : double black_hole_mass_{std::numeric_limits<double>::signaling_NaN()};
164 0 : double black_hole_spin_{std::numeric_limits<double>::signaling_NaN()};
165 0 : double orbital_radius_{std::numeric_limits<double>::signaling_NaN()};
166 0 : int m_mode_number_{};
167 0 : std::optional<std::array<double, 4>> hyperboloidal_slicing_transitions_{};
168 0 : bool penetrating_horizon_{false};
169 : };
170 :
171 0 : bool operator!=(const CircularOrbit& lhs, const CircularOrbit& rhs);
172 :
173 : } // namespace GrSelfForce::AnalyticData
|
