Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cmath>
7 : #include <cstddef>
8 :
9 : #include "DataStructures/DataBox/Tag.hpp"
10 : #include "DataStructures/DataVector.hpp"
11 : #include "DataStructures/TaggedTuple.hpp"
12 : #include "DataStructures/Tensor/TypeAliases.hpp"
13 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
14 : #include "NumericalAlgorithms/Spectral/CollocationPoints.hpp"
15 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
16 : #include "Options/Context.hpp"
17 : #include "Options/String.hpp"
18 : #include "ParallelAlgorithms/Interpolation/Protocols/ComputeTargetPoints.hpp"
19 : #include "ParallelAlgorithms/Interpolation/Tags.hpp"
20 : #include "Utilities/PrettyType.hpp"
21 : #include "Utilities/ProtocolHelpers.hpp"
22 : #include "Utilities/Requires.hpp"
23 : #include "Utilities/TMPL.hpp"
24 :
25 : /// \cond
26 : namespace PUP {
27 : class er;
28 : } // namespace PUP
29 : namespace db {
30 : template <typename TagsList>
31 : class DataBox;
32 : } // namespace db
33 : namespace intrp {
34 : namespace Tags {
35 : template <typename TemporalId>
36 : struct TemporalIds;
37 : } // namespace Tags
38 : } // namespace intrp
39 : /// \endcond
40 :
41 : namespace intrp {
42 :
43 : namespace OptionHolders {
44 : /// \brief A solid torus of points, useful, e.g., when measuring data from an
45 : /// accretion disc.
46 : ///
47 : /// The torus's cross section (e.g., a cut at \f$\phi=0\f$) is a wedge-like
48 : /// shape bounded by \f$r_{\text{min}} \le r \le r_{\text{max}}\f$ and
49 : /// \f$\theta_{\text{min}} \le \theta \le \theta_{\text{max}}\f$.
50 : ///
51 : /// The grid points are located on surfaces of constant \f$r\f$, \f$\theta\f$,
52 : /// and \f$\phi\f$. There are `NumberRadialPoints` points in the radial
53 : /// direction between `MinRadius` and `MaxRadius` (including these endpoints);
54 : /// `NumberThetaPoints` points in the \f$\theta\f$ direction between `MinTheta`
55 : /// and `MaxTheta` (including these endpoints); `NumberPhiPoints` points in the
56 : /// \f$\phi\f$ direction (with one point always at \f$\phi=0\f$).
57 : ///
58 : /// By default, the points follow a Legendre Gauss-Lobatto distribution in the
59 : /// \f$r\f$ and \f$\theta\f$ directions, and a uniform distribution in the
60 : /// \f$\phi\f$ direction. The distribution in the \f$r\f$ (and/or \f$\theta\f$)
61 : /// direction can be made uniform using the `UniformRadialGrid` (and/or
62 : /// `UniformThetaGrid`) option.
63 : ///
64 : /// The `target_points` form a 3D mesh ordered with \f$r\f$ varying fastest,
65 : /// then \f$\theta\f$, and finally \f$\phi\f$ varying slowest.
66 : ///
67 : /// \note Input coordinates (radii, angles) are interpreted in the
68 : /// `Frame::Inertial`
69 1 : struct WedgeSectionTorus {
70 0 : struct MinRadius {
71 0 : using type = double;
72 0 : static constexpr Options::String help = {"Inner radius of torus"};
73 0 : static type lower_bound() { return 0.0; }
74 : };
75 0 : struct MaxRadius {
76 0 : using type = double;
77 0 : static constexpr Options::String help = {"Outer radius of torus"};
78 0 : static type lower_bound() { return 0.0; }
79 : };
80 0 : struct MinTheta {
81 0 : using type = double;
82 0 : static constexpr Options::String help = {"Angle of top of wedge (radians)"};
83 0 : static type lower_bound() { return 0.0; }
84 0 : static type upper_bound() { return M_PI; }
85 : };
86 0 : struct MaxTheta {
87 0 : using type = double;
88 0 : static constexpr Options::String help = {
89 : "Angle of bottom of wedge (radians)"};
90 0 : static type lower_bound() { return 0.0; }
91 0 : static type upper_bound() { return M_PI; }
92 : };
93 0 : struct NumberRadialPoints {
94 0 : using type = size_t;
95 0 : static constexpr Options::String help = {
96 : "Number of radial points, including endpoints"};
97 0 : static type lower_bound() { return 2; }
98 : };
99 0 : struct NumberThetaPoints {
100 0 : using type = size_t;
101 0 : static constexpr Options::String help = {
102 : "Number of theta points, including endpoints"};
103 0 : static type lower_bound() { return 2; }
104 : };
105 0 : struct NumberPhiPoints {
106 0 : using type = size_t;
107 0 : static constexpr Options::String help = {"Number of phi points"};
108 0 : static type lower_bound() { return 1; }
109 : };
110 0 : struct UniformRadialGrid {
111 0 : using type = bool;
112 0 : static constexpr Options::String help = {"Use uniform radial grid"};
113 : };
114 0 : struct UniformThetaGrid {
115 0 : using type = bool;
116 0 : static constexpr Options::String help = {"Use uniform theta grid"};
117 : };
118 :
119 0 : using options =
120 : tmpl::list<MinRadius, MaxRadius, MinTheta, MaxTheta, NumberRadialPoints,
121 : NumberThetaPoints, NumberPhiPoints, UniformRadialGrid,
122 : UniformThetaGrid>;
123 0 : static constexpr Options::String help = {
124 : "A torus extending from MinRadius to MaxRadius in r, MinTheta to MaxTheta"
125 : " in theta, and 2pi in phi."};
126 :
127 0 : WedgeSectionTorus(double min_radius_in, double max_radius_in,
128 : double min_theta_in, double max_theta_in,
129 : size_t number_of_radial_points_in,
130 : size_t number_of_theta_points_in,
131 : size_t number_of_phi_points_in,
132 : bool use_uniform_radial_grid_in,
133 : bool use_uniform_theta_grid_in,
134 : const Options::Context& context = {});
135 :
136 0 : WedgeSectionTorus() = default;
137 0 : WedgeSectionTorus(const WedgeSectionTorus& /*rhs*/) = default;
138 0 : WedgeSectionTorus& operator=(const WedgeSectionTorus& /*rhs*/) = delete;
139 0 : WedgeSectionTorus(WedgeSectionTorus&& /*rhs*/) = default;
140 0 : WedgeSectionTorus& operator=(WedgeSectionTorus&& /*rhs*/) = default;
141 0 : ~WedgeSectionTorus() = default;
142 :
143 : // NOLINTNEXTLINE(google-runtime-references)
144 0 : void pup(PUP::er& p);
145 :
146 0 : double min_radius;
147 0 : double max_radius;
148 0 : double min_theta;
149 0 : double max_theta;
150 0 : size_t number_of_radial_points;
151 0 : size_t number_of_theta_points;
152 0 : size_t number_of_phi_points;
153 0 : bool use_uniform_radial_grid;
154 0 : bool use_uniform_theta_grid;
155 : };
156 :
157 0 : bool operator==(const WedgeSectionTorus& lhs, const WedgeSectionTorus& rhs);
158 0 : bool operator!=(const WedgeSectionTorus& lhs, const WedgeSectionTorus& rhs);
159 :
160 : } // namespace OptionHolders
161 :
162 : namespace OptionTags {
163 : template <typename InterpolationTargetTag>
164 0 : struct WedgeSectionTorus {
165 0 : using type = OptionHolders::WedgeSectionTorus;
166 0 : static constexpr Options::String help{
167 : "Options for interpolation onto Kerr horizon."};
168 0 : static std::string name() {
169 : return pretty_type::name<InterpolationTargetTag>();
170 : }
171 0 : using group = InterpolationTargets;
172 : };
173 : } // namespace OptionTags
174 :
175 : namespace Tags {
176 : template <typename InterpolationTargetTag>
177 0 : struct WedgeSectionTorus : db::SimpleTag {
178 0 : using type = OptionHolders::WedgeSectionTorus;
179 0 : using option_tags =
180 : tmpl::list<OptionTags::WedgeSectionTorus<InterpolationTargetTag>>;
181 :
182 0 : static constexpr bool pass_metavariables = false;
183 0 : static type create_from_options(const type& option) { return option; }
184 : };
185 : } // namespace Tags
186 :
187 : namespace TargetPoints {
188 : /// \brief Computes points in a wedge-sectioned torus.
189 : ///
190 : /// Conforms to the intrp::protocols::ComputeTargetPoints protocol
191 : ///
192 : /// For requirements on InterpolationTargetTag, see
193 : /// intrp::protocols::InterpolationTargetTag
194 : template <typename InterpolationTargetTag>
195 1 : struct WedgeSectionTorus
196 : : tt::ConformsTo<intrp::protocols::ComputeTargetPoints> {
197 0 : using const_global_cache_tags =
198 : tmpl::list<Tags::WedgeSectionTorus<InterpolationTargetTag>>;
199 0 : using is_sequential = std::false_type;
200 0 : using frame = Frame::Inertial;
201 :
202 : template <typename Metavariables, typename DbTags>
203 0 : static tnsr::I<DataVector, 3, Frame::Inertial> points(
204 : const db::DataBox<DbTags>& box,
205 : const tmpl::type_<Metavariables>& /*meta*/) {
206 : const auto& options =
207 : get<Tags::WedgeSectionTorus<InterpolationTargetTag>>(box);
208 :
209 : // Compute locations of constant r/theta/phi surfaces
210 : const size_t num_radial = options.number_of_radial_points;
211 : const DataVector radii_1d = [&num_radial, &options]() {
212 : DataVector result(num_radial);
213 : if (options.use_uniform_radial_grid) {
214 : // uniform point distribution
215 : const double coefficient =
216 : (options.max_radius - options.min_radius) / (num_radial - 1.0);
217 : for (size_t r = 0; r < num_radial; ++r) {
218 : result[r] = options.min_radius + coefficient * r;
219 : }
220 : } else {
221 : // radial Legendre-Gauss-Lobatto distribution via linear rescaling
222 : const double mean = 0.5 * (options.max_radius + options.min_radius);
223 : const double diff = 0.5 * (options.max_radius - options.min_radius);
224 : result =
225 : mean + diff * Spectral::collocation_points<
226 : Spectral::Basis::Legendre,
227 : Spectral::Quadrature::GaussLobatto>(num_radial);
228 : }
229 : return result;
230 : }();
231 : const size_t num_theta = options.number_of_theta_points;
232 : const DataVector thetas_1d = [&num_theta, &options]() {
233 : DataVector result(num_theta);
234 : if (options.use_uniform_theta_grid) {
235 : // uniform point distribution
236 : const double coefficient =
237 : (options.max_theta - options.min_theta) / (num_theta - 1.0);
238 : for (size_t theta = 0; theta < num_theta; ++theta) {
239 : result[theta] = options.min_theta + coefficient * theta;
240 : }
241 : } else {
242 : // Legendre-Gauss-Lobatto point distribution (in theta)
243 : const double mean = 0.5 * (options.max_theta + options.min_theta);
244 : const double diff = 0.5 * (options.max_theta - options.min_theta);
245 : result =
246 : mean + diff * Spectral::collocation_points<
247 : Spectral::Basis::Legendre,
248 : Spectral::Quadrature::GaussLobatto>(num_theta);
249 : }
250 : return result;
251 : }();
252 : const size_t num_phi = options.number_of_phi_points;
253 : const DataVector phis_1d = [&num_phi]() {
254 : DataVector result(num_phi);
255 : // We do NOT want a grid point at phi = 2pi, as this would duplicate the
256 : // phi = 0 data. So, divide by num_phi rather than (n-1) as elsewhere.
257 : const double coefficient = 2.0 * M_PI / num_phi;
258 : for (size_t phi = 0; phi < num_phi; ++phi) {
259 : result[phi] = coefficient * phi;
260 : }
261 : return result;
262 : }();
263 :
264 : // Take tensor product to get full 3D r/theta/phi points
265 : const size_t num_total = num_radial * num_theta * num_phi;
266 : DataVector radii(num_total);
267 : DataVector thetas(num_total);
268 : DataVector phis(num_total);
269 : for (size_t phi = 0; phi < num_phi; ++phi) {
270 : for (size_t theta = 0; theta < num_theta; ++theta) {
271 : for (size_t r = 0; r < num_radial; ++r) {
272 : const size_t i =
273 : r + theta * num_radial + phi * num_theta * num_radial;
274 : radii[i] = radii_1d[r];
275 : thetas[i] = thetas_1d[theta];
276 : phis[i] = phis_1d[phi];
277 : }
278 : }
279 : }
280 :
281 : // Compute x/y/z coordinates
282 : // Note: theta measured from +z axis, phi measured from +x axis
283 : tnsr::I<DataVector, 3, Frame::Inertial> target_points(num_total);
284 : get<0>(target_points) = radii * sin(thetas) * cos(phis);
285 : get<1>(target_points) = radii * sin(thetas) * sin(phis);
286 : get<2>(target_points) = radii * cos(thetas);
287 :
288 : return target_points;
289 : }
290 :
291 : template <typename Metavariables, typename DbTags, typename TemporalId>
292 0 : static tnsr::I<DataVector, 3, Frame::Inertial> points(
293 : const db::DataBox<DbTags>& box, const tmpl::type_<Metavariables>& meta,
294 : const TemporalId& /*temporal_id*/) {
295 : return points(box, meta);
296 : }
297 : };
298 :
299 : } // namespace TargetPoints
300 : } // namespace intrp
|
