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 :
9 : #include "DataStructures/DataBox/DataBox.hpp"
10 : #include "DataStructures/DataBox/Tag.hpp"
11 : #include "DataStructures/TaggedTuple.hpp"
12 : #include "DataStructures/Tensor/TypeAliases.hpp"
13 : #include "DataStructures/Transpose.hpp"
14 : #include "NumericalAlgorithms/SphericalHarmonics/AngularOrdering.hpp"
15 : #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
16 : #include "NumericalAlgorithms/Strahlkorper/Strahlkorper.hpp"
17 : #include "NumericalAlgorithms/Strahlkorper/Tags.hpp"
18 : #include "Options/String.hpp"
19 : #include "Parallel/GlobalCache.hpp"
20 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
21 : #include "ParallelAlgorithms/Interpolation/Protocols/ComputeTargetPoints.hpp"
22 : #include "ParallelAlgorithms/Interpolation/Tags.hpp"
23 : #include "PointwiseFunctions/GeneralRelativity/KerrHorizon.hpp"
24 : #include "Utilities/PrettyType.hpp"
25 : #include "Utilities/ProtocolHelpers.hpp"
26 : #include "Utilities/Requires.hpp"
27 : #include "Utilities/TMPL.hpp"
28 :
29 : /// \cond
30 : class DataVector;
31 : namespace PUP {
32 : class er;
33 : } // namespace PUP
34 : namespace db {
35 : template <typename TagsList>
36 : class DataBox;
37 : } // namespace db
38 : namespace intrp {
39 : namespace Tags {
40 : template <typename TemporalId>
41 : struct TemporalIds;
42 : } // namespace Tags
43 : } // namespace intrp
44 : /// \endcond
45 :
46 : namespace intrp {
47 :
48 0 : namespace OptionHolders {
49 : /// A surface that conforms to the horizon of a Kerr black hole in
50 : /// Kerr-Schild coordinates.
51 : ///
52 : /// \details In addition to the parameters for the Kerr black hole, this holder
53 : /// contains the `LMax` which encodes the angular resolution of the spherical
54 : /// harmonic basis and `AngularOrdering` which encodes the collocation
55 : /// ordering.
56 1 : struct KerrHorizon {
57 0 : struct LMax {
58 0 : using type = size_t;
59 0 : static constexpr Options::String help = {
60 : "KerrHorizon is expanded in Ylms up to l=LMax"};
61 : };
62 0 : struct Center {
63 0 : using type = std::array<double, 3>;
64 0 : static constexpr Options::String help = {"Center of black hole"};
65 : };
66 0 : struct Mass {
67 0 : using type = double;
68 0 : static constexpr Options::String help = {"Mass of black hole"};
69 : };
70 0 : struct DimensionlessSpin {
71 0 : using type = std::array<double, 3>;
72 0 : static constexpr Options::String help = {
73 : "Dimensionless spin of black hole"};
74 : };
75 0 : struct AngularOrdering {
76 0 : using type = ylm::AngularOrdering;
77 0 : static constexpr Options::String help = {
78 : "Chooses theta,phi ordering in 2d array"};
79 : };
80 0 : using options =
81 : tmpl::list<LMax, Center, Mass, DimensionlessSpin, AngularOrdering>;
82 0 : static constexpr Options::String help = {
83 : "A Strahlkorper conforming to the horizon (in Kerr-Schild coordinates)"
84 : " of a Kerr black hole with a specified center, mass, and spin."};
85 :
86 0 : KerrHorizon(size_t l_max_in, std::array<double, 3> center_in, double mass_in,
87 : std::array<double, 3> dimensionless_spin_in,
88 : ylm::AngularOrdering angular_ordering_in,
89 : const Options::Context& context = {});
90 :
91 0 : KerrHorizon() = default;
92 0 : KerrHorizon(const KerrHorizon& /*rhs*/) = default;
93 0 : KerrHorizon& operator=(const KerrHorizon& /*rhs*/) = delete;
94 0 : KerrHorizon(KerrHorizon&& /*rhs*/) = default;
95 0 : KerrHorizon& operator=(KerrHorizon&& /*rhs*/) = default;
96 0 : ~KerrHorizon() = default;
97 :
98 : // NOLINTNEXTLINE(google-runtime-references)
99 0 : void pup(PUP::er& p);
100 :
101 0 : size_t l_max{};
102 0 : std::array<double, 3> center{};
103 0 : double mass{};
104 0 : std::array<double, 3> dimensionless_spin{};
105 0 : ylm::AngularOrdering angular_ordering;
106 : };
107 :
108 0 : bool operator==(const KerrHorizon& lhs, const KerrHorizon& rhs);
109 0 : bool operator!=(const KerrHorizon& lhs, const KerrHorizon& rhs);
110 :
111 : } // namespace OptionHolders
112 :
113 : namespace OptionTags {
114 : template <typename InterpolationTargetTag>
115 0 : struct KerrHorizon {
116 0 : using type = OptionHolders::KerrHorizon;
117 0 : static constexpr Options::String help{
118 : "Options for interpolation onto Kerr horizon."};
119 0 : static std::string name() {
120 : return pretty_type::name<InterpolationTargetTag>();
121 : }
122 0 : using group = InterpolationTargets;
123 : };
124 : } // namespace OptionTags
125 :
126 : namespace Tags {
127 : template <typename InterpolationTargetTag>
128 0 : struct KerrHorizon : db::SimpleTag {
129 0 : using type = OptionHolders::KerrHorizon;
130 0 : using option_tags =
131 : tmpl::list<OptionTags::KerrHorizon<InterpolationTargetTag>>;
132 :
133 0 : static constexpr bool pass_metavariables = false;
134 0 : static type create_from_options(const type& option) { return option; }
135 : };
136 : } // namespace Tags
137 :
138 0 : namespace TargetPoints {
139 : /// \brief Computes points on a Kerr horizon.
140 : ///
141 : /// The points are such that they conform to the horizon of a Kerr
142 : /// black hole (in Kerr-Schild coordinates) with given center, mass,
143 : /// and dimensionless spin, as specified in the options.
144 : ///
145 : /// \note The returned points are not actually on the horizon;
146 : /// instead, they are collocation points of a Strahlkorper whose
147 : /// spectral representation matches the horizon shape up to order
148 : /// LMax, where LMax is the spherical-harmonic order specified in the
149 : /// options. As LMax increases, the returned points converge to the
150 : /// horizon.
151 : ///
152 : /// Conforms to the intrp::protocols::ComputeTargetPoints protocol
153 : ///
154 : /// For requirements on InterpolationTargetTag, see
155 : /// intrp::protocols::InterpolationTargetTag
156 : template <typename InterpolationTargetTag, typename Frame>
157 1 : struct KerrHorizon : tt::ConformsTo<intrp::protocols::ComputeTargetPoints> {
158 0 : using const_global_cache_tags =
159 : tmpl::list<Tags::KerrHorizon<InterpolationTargetTag>>;
160 0 : using is_sequential = std::false_type;
161 0 : using frame = Frame;
162 :
163 0 : using simple_tags = typename ylm::Tags::items_tags<Frame>;
164 0 : using compute_tags = typename ylm::Tags::compute_items_tags<Frame>;
165 :
166 : template <typename DbTags, typename Metavariables>
167 0 : static void initialize(const gsl::not_null<db::DataBox<DbTags>*> box,
168 : const Parallel::GlobalCache<Metavariables>& cache) {
169 : const auto& kerr_horizon =
170 : Parallel::get<Tags::KerrHorizon<InterpolationTargetTag>>(cache);
171 :
172 : // Make a Strahlkorper with the correct shape.
173 : ylm::Strahlkorper<Frame> strahlkorper(
174 : kerr_horizon.l_max, kerr_horizon.l_max,
175 : get(gr::Solutions::kerr_horizon_radius(
176 : ::ylm::Spherepack(kerr_horizon.l_max, kerr_horizon.l_max)
177 : .theta_phi_points(),
178 : kerr_horizon.mass, kerr_horizon.dimensionless_spin)),
179 : kerr_horizon.center);
180 : Initialization::mutate_assign<simple_tags>(box, std::move(strahlkorper));
181 : }
182 :
183 : template <typename Metavariables, typename DbTags, typename TemporalId>
184 0 : static tnsr::I<DataVector, 3, Frame> points(
185 : const db::DataBox<DbTags>& box,
186 : const tmpl::type_<Metavariables>& metavariables_v,
187 : const TemporalId& /*temporal_id*/) {
188 : return points(box, metavariables_v);
189 : }
190 : template <typename Metavariables, typename DbTags>
191 0 : static tnsr::I<DataVector, 3, Frame> points(
192 : const db::DataBox<DbTags>& box,
193 : const tmpl::type_<Metavariables>& /*meta*/) {
194 : const auto& kerr_horizon =
195 : db::get<Tags::KerrHorizon<InterpolationTargetTag>>(box);
196 : if (kerr_horizon.angular_ordering == ylm::AngularOrdering::Strahlkorper) {
197 : return db::get<ylm::Tags::CartesianCoords<Frame>>(box);
198 : } else {
199 : const auto& strahlkorper = db::get<ylm::Tags::Strahlkorper<Frame>>(box);
200 : const auto& coords = db::get<ylm::Tags::CartesianCoords<Frame>>(box);
201 : const auto physical_extents =
202 : strahlkorper.ylm_spherepack().physical_extents();
203 : auto transposed_coords =
204 : tnsr::I<DataVector, 3, Frame>(get<0>(coords).size());
205 : for (size_t i = 0; i < 3; ++i) {
206 : transpose(make_not_null(&transposed_coords.get(i)), coords.get(i),
207 : physical_extents[0], physical_extents[1]);
208 : }
209 : return transposed_coords;
210 : }
211 : }
212 : };
213 :
214 : } // namespace TargetPoints
215 : } // namespace intrp
|
