Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <array>
8 : #include <cstddef>
9 : #include <limits>
10 : #include <set>
11 : #include <variant>
12 : #include <vector>
13 :
14 : #include "DataStructures/DataBox/DataBox.hpp"
15 : #include "DataStructures/Transpose.hpp"
16 : #include "NumericalAlgorithms/SphericalHarmonics/AngularOrdering.hpp"
17 : #include "NumericalAlgorithms/Strahlkorper/Strahlkorper.hpp"
18 : #include "NumericalAlgorithms/Strahlkorper/Tags.hpp"
19 : #include "Options/String.hpp"
20 : #include "Parallel/GlobalCache.hpp"
21 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
22 : #include "ParallelAlgorithms/Interpolation/Protocols/ComputeTargetPoints.hpp"
23 : #include "ParallelAlgorithms/Interpolation/Tags.hpp"
24 : #include "Utilities/Algorithm.hpp"
25 : #include "Utilities/ErrorHandling/Assert.hpp"
26 : #include "Utilities/PrettyType.hpp"
27 : #include "Utilities/ProtocolHelpers.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 : namespace OptionHolders {
49 : /// A series of concentric spherical surfaces.
50 : ///
51 : /// \details The parameter `LMax` sets the number of collocation points on
52 : /// each spherical surface equal to `(l_max + 1) * (2 * l_max + 1)`. The
53 : /// parameter `AngularOrdering` encodes the collocation ordering. For example,
54 : /// the apparent horizon finder relies on spherepack routines that require
55 : /// `Strahlkorper` for `AngularOrdering`, and using these surfaces for a CCE
56 : /// worldtube requires `Cce` for `AngularOrdering`.
57 1 : struct Sphere {
58 0 : struct LMax {
59 0 : using type = size_t;
60 0 : static constexpr Options::String help = {
61 : "The number of collocation points on each sphere will be equal to "
62 : "`(l_max + 1) * (2 * l_max + 1)`"};
63 : };
64 0 : struct Center {
65 0 : using type = std::array<double, 3>;
66 0 : static constexpr Options::String help = {"Center of every sphere"};
67 : };
68 0 : struct Radius {
69 0 : using type = std::variant<double, std::vector<double>>;
70 0 : static constexpr Options::String help = {"Radius of the sphere(s)"};
71 : };
72 0 : struct AngularOrdering {
73 0 : using type = ylm::AngularOrdering;
74 0 : static constexpr Options::String help = {
75 : "Chooses theta,phi ordering in 2d array"};
76 : };
77 0 : using options = tmpl::list<LMax, Center, Radius, AngularOrdering>;
78 0 : static constexpr Options::String help = {
79 : "An arbitrary number of spherical surface."};
80 0 : Sphere(const size_t l_max_in, const std::array<double, 3> center_in,
81 : const typename Radius::type& radius_in,
82 : ylm::AngularOrdering angular_ordering_in,
83 : const Options::Context& context = {});
84 :
85 0 : Sphere() = default;
86 :
87 : // NOLINTNEXTLINE(google-runtime-references)
88 0 : void pup(PUP::er& p);
89 :
90 0 : size_t l_max{0};
91 0 : std::array<double, 3> center{std::numeric_limits<double>::signaling_NaN()};
92 0 : std::set<double> radii;
93 0 : ylm::AngularOrdering angular_ordering;
94 : };
95 :
96 0 : bool operator==(const Sphere& lhs, const Sphere& rhs);
97 0 : bool operator!=(const Sphere& lhs, const Sphere& rhs);
98 :
99 : } // namespace OptionHolders
100 :
101 : namespace OptionTags {
102 : template <typename InterpolationTargetTag>
103 0 : struct Sphere {
104 0 : using type = OptionHolders::Sphere;
105 0 : static constexpr Options::String help{
106 : "Options for interpolation onto a sphere(s)."};
107 0 : static std::string name() {
108 : return pretty_type::name<InterpolationTargetTag>();
109 : }
110 0 : using group = InterpolationTargets;
111 : };
112 : } // namespace OptionTags
113 :
114 : namespace Tags {
115 : template <typename InterpolationTargetTag>
116 0 : struct Sphere : db::SimpleTag {
117 0 : using type = OptionHolders::Sphere;
118 0 : using option_tags = tmpl::list<OptionTags::Sphere<InterpolationTargetTag>>;
119 :
120 0 : static constexpr bool pass_metavariables = false;
121 0 : static OptionHolders::Sphere create_from_options(
122 : const OptionHolders::Sphere& option) {
123 : return option;
124 : }
125 : };
126 :
127 : template <typename Frame>
128 0 : struct AllCoords : db::SimpleTag {
129 0 : using type = tnsr::I<DataVector, 3, Frame>;
130 : };
131 : } // namespace Tags
132 :
133 : namespace TargetPoints {
134 : /// \brief Computes points on spherical surfaces.
135 : ///
136 : /// Conforms to the intrp::protocols::ComputeTargetPoints protocol
137 : ///
138 : /// For requirements on InterpolationTargetTag, see
139 : /// intrp::protocols::InterpolationTargetTag
140 : template <typename InterpolationTargetTag, typename Frame>
141 1 : struct Sphere : tt::ConformsTo<intrp::protocols::ComputeTargetPoints> {
142 0 : using const_global_cache_tags =
143 : tmpl::list<Tags::Sphere<InterpolationTargetTag>>;
144 0 : using is_sequential = std::false_type;
145 0 : using frame = Frame;
146 :
147 0 : using simple_tags =
148 : tmpl::list<ylm::Tags::Strahlkorper<Frame>, Tags::AllCoords<Frame>>;
149 0 : using compute_tags = typename ylm::Tags::compute_items_tags<Frame>;
150 :
151 : template <typename DbTags, typename Metavariables>
152 0 : static void initialize(const gsl::not_null<db::DataBox<DbTags>*> box,
153 : const Parallel::GlobalCache<Metavariables>& cache) {
154 : const auto& sphere =
155 : Parallel::get<Tags::Sphere<InterpolationTargetTag>>(cache);
156 : const size_t l_max = sphere.l_max;
157 : const auto& radii = sphere.radii;
158 :
159 : // Total number of points is number of points for one sphere times the
160 : // number of spheres we use.
161 : const size_t num_points = radii.size() * (l_max + 1) * (2 * l_max + 1);
162 :
163 : tnsr::I<DataVector, 3, Frame> all_coords{num_points};
164 :
165 : size_t index = 0;
166 : for (const double radius : radii) {
167 : ylm::Strahlkorper<Frame> strahlkorper(
168 : l_max, l_max, DataVector{(l_max + 1) * (2 * l_max + 1), radius},
169 : sphere.center);
170 :
171 : db::mutate<ylm::Tags::Strahlkorper<Frame>>(
172 : [&strahlkorper](const gsl::not_null<ylm::Strahlkorper<Frame>*>
173 : local_strahlkorper) {
174 : *local_strahlkorper = std::move(strahlkorper);
175 : },
176 : box);
177 :
178 : // This copy is ok because it's just in initialization
179 : auto coords = db::get<ylm::Tags::CartesianCoords<Frame>>(*box);
180 :
181 : // If the angular ordering is Strahlkorper then we don't have to do
182 : // anything to the coords because they are already in the right order
183 : if (sphere.angular_ordering == ylm::AngularOrdering::Cce) {
184 : const auto physical_extents =
185 : strahlkorper.ylm_spherepack().physical_extents();
186 : auto transposed_coords =
187 : tnsr::I<DataVector, 3, Frame>(get<0>(coords).size());
188 :
189 : for (size_t i = 0; i < 3; ++i) {
190 : transpose(make_not_null(&transposed_coords.get(i)), coords.get(i),
191 : physical_extents[0], physical_extents[1]);
192 : }
193 : coords = std::move(transposed_coords);
194 : }
195 :
196 : const size_t tmp_index = index;
197 : for (size_t i = 0; i < 3; ++i) {
198 : for (size_t local_index = 0; local_index < coords.get(i).size();
199 : local_index++, index++) {
200 : all_coords.get(i)[index] = coords.get(i)[local_index];
201 : }
202 : if (i != 2) {
203 : index = tmp_index;
204 : }
205 : }
206 : }
207 :
208 : // If this fails, there is a bug. Can't really test it
209 : // LCOV_EXCL_START
210 : ASSERT(index == all_coords.get(0).size(),
211 : "Didn't initialize points of Sphere target correctly. index = "
212 : << index << " all_coords.size() = " << all_coords.get(0).size());
213 : // LCOV_EXCL_STOP
214 :
215 : Initialization::mutate_assign<tmpl::list<Tags::AllCoords<Frame>>>(
216 : box, std::move(all_coords));
217 : }
218 :
219 : template <typename Metavariables, typename DbTags, typename TemporalId>
220 0 : static tnsr::I<DataVector, 3, Frame> points(
221 : const db::DataBox<DbTags>& box,
222 : const tmpl::type_<Metavariables>& metavariables_v,
223 : const TemporalId& /*temporal_id*/) {
224 : return points(box, metavariables_v);
225 : }
226 : template <typename Metavariables, typename DbTags>
227 0 : static tnsr::I<DataVector, 3, Frame> points(
228 : const db::DataBox<DbTags>& box,
229 : const tmpl::type_<Metavariables>& /*meta*/) {
230 : return db::get<Tags::AllCoords<Frame>>(box);
231 : }
232 : };
233 : } // namespace TargetPoints
234 : } // namespace intrp
|
