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 <limits>
8 : #include <memory>
9 : #include <string>
10 : #include <type_traits>
11 : #include <unordered_map>
12 :
13 : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
14 : #include "Domain/CoordinateMaps/Identity.hpp"
15 : #include "Domain/CoordinateMaps/TimeDependent/CubicScale.hpp"
16 : #include "Domain/CoordinateMaps/TimeDependent/RotScaleTrans.hpp"
17 : #include "Domain/CoordinateMaps/TimeDependent/Rotation.hpp"
18 : #include "Domain/CoordinateMaps/TimeDependent/Shape.hpp"
19 : #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp"
20 : #include "Domain/Creators/ShapeMapOptions.hpp"
21 : #include "Domain/Creators/SphereTimeDependentMaps.hpp"
22 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
23 : #include "Domain/Structure/ObjectLabel.hpp"
24 : #include "Options/Auto.hpp"
25 : #include "Options/Context.hpp"
26 : #include "Options/Options.hpp"
27 : #include "Options/String.hpp"
28 : #include "Utilities/GetOutput.hpp"
29 : #include "Utilities/TMPL.hpp"
30 :
31 : /// \cond
32 : namespace Frame {
33 : struct Grid;
34 : struct Distorted;
35 : struct Inertial;
36 : } // namespace Frame
37 : /// \endcond
38 :
39 : /// Namespace used to hold things used in both the BinaryCompactObject and
40 : /// CylindricalBinaryCompactObject domain creators.
41 1 : namespace domain::creators::bco {
42 :
43 : /*!
44 : * \brief Create a set of centers of objects for the binary domains.
45 : *
46 : * \details Will add the following centers to the set:
47 : *
48 : * - Center: The origin
49 : * - CenterA: Center of object A
50 : * - CenterB: Center of object B
51 : *
52 : * \return Object required by the DomainCreator%s
53 : */
54 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
55 1 : create_grid_anchors(const std::array<double, 3>& center_a,
56 : const std::array<double, 3>& center_b);
57 :
58 : namespace detail {
59 : // Convenience type alias
60 : template <typename... Maps>
61 : using gi_map = domain::CoordinateMap<Frame::Grid, Frame::Inertial, Maps...>;
62 : template <typename... Maps>
63 : using gd_map = domain::CoordinateMap<Frame::Grid, Frame::Distorted, Maps...>;
64 : template <typename... Maps>
65 : using di_map =
66 : domain::CoordinateMap<Frame::Distorted, Frame::Inertial, Maps...>;
67 :
68 : template <typename List>
69 : struct power_set {
70 : using rest = typename power_set<tmpl::pop_front<List>>::type;
71 : using type = tmpl::append<
72 : rest, tmpl::transform<rest, tmpl::lazy::push_front<
73 : tmpl::_1, tmpl::pin<tmpl::front<List>>>>>;
74 : };
75 :
76 : template <>
77 : struct power_set<tmpl::list<>> {
78 : using type = tmpl::list<tmpl::list<>>;
79 : };
80 :
81 : template <typename SourceFrame, typename TargetFrame, typename Maps>
82 : using produce_all_maps_helper =
83 : tmpl::wrap<tmpl::push_front<Maps, SourceFrame, TargetFrame>,
84 : domain::CoordinateMap>;
85 :
86 : /*
87 : * This will produce all the possible combinations of maps we will need for the
88 : * maps_list. It does so using a power set. Say you have 3 maps, Map1, Map2,
89 : * Map3. The end result of this will be a list with the following map
90 : * combinations where these combinations may not appear in the same position in
91 : * the list, but the order of each combination is fixed:
92 : * - Map1
93 : * - Map2
94 : * - Map3
95 : * - Map1, Map2
96 : * - Map1, Map3
97 : * - Map2, Map3
98 : * - Map1, Map2, Map3
99 : */
100 : template <typename SourceFrame, typename TargetFrame, typename... Maps>
101 : using produce_all_maps = tmpl::transform<
102 : tmpl::remove<typename power_set<tmpl::list<Maps...>>::type, tmpl::list<>>,
103 : tmpl::bind<produce_all_maps_helper, tmpl::pin<SourceFrame>,
104 : tmpl::pin<TargetFrame>, tmpl::_1>>;
105 : } // namespace detail
106 :
107 : /*!
108 : * \brief This holds all options related to the time dependent maps of the
109 : * binary compact object domains.
110 : *
111 : * \details Since both domains will have the same (overall) time dependent maps,
112 : * their options are going to be the same as well. The options won't be exactly
113 : * the same though, so there is a \p IsCylindrical template parameter to
114 : * distinguish.
115 : *
116 : * This class will create the FunctionsOfTime needed for the binary compact
117 : * object domains as well as the actual CoordinateMap%s themselves
118 : *
119 : * \note This struct contains no information about what blocks the time
120 : * dependent maps will go in.
121 : */
122 : template <bool IsCylindrical>
123 1 : struct TimeDependentMapOptions {
124 : private:
125 : template <typename SourceFrame, typename TargetFrame>
126 0 : using MapType =
127 : std::unique_ptr<domain::CoordinateMapBase<SourceFrame, TargetFrame, 3>>;
128 : // Time-dependent maps
129 0 : using Expansion = domain::CoordinateMaps::TimeDependent::CubicScale<3>;
130 0 : using Rotation = domain::CoordinateMaps::TimeDependent::Rotation<3>;
131 0 : using RotScaleTrans = domain::CoordinateMaps::TimeDependent::RotScaleTrans<3>;
132 0 : using Shape = domain::CoordinateMaps::TimeDependent::Shape;
133 0 : using Identity = domain::CoordinateMaps::Identity<3>;
134 :
135 : public:
136 0 : using maps_list = tmpl::append<
137 : // We need this odd-one-out Identity map because all maps are optional to
138 : // specify. It's possible to specify a shape map, but no expansion or
139 : // rotation map. So we have a grid to distorted map, and need an identity
140 : // distorted to inertial map. We don't need an identity grid to distorted
141 : // map because if a user requests the grid to distorted frame map with a
142 : // distorted frame, but didn't specify shape map options, an error occurs.
143 : tmpl::list<detail::di_map<Identity>>,
144 : detail::produce_all_maps<Frame::Grid, Frame::Inertial, Shape, Expansion,
145 : Rotation>,
146 : detail::produce_all_maps<Frame::Grid, Frame::Distorted, Shape>,
147 : detail::produce_all_maps<Frame::Distorted, Frame::Inertial, Expansion,
148 : Rotation>,
149 : detail::produce_all_maps<Frame::Grid, Frame::Inertial, Shape,
150 : RotScaleTrans>,
151 : detail::produce_all_maps<Frame::Distorted, Frame::Inertial,
152 : RotScaleTrans>>;
153 :
154 : /// \brief The initial time of the functions of time.
155 1 : struct InitialTime {
156 0 : using type = double;
157 0 : static constexpr Options::String help = {
158 : "The initial time of the functions of time"};
159 : };
160 :
161 : /// \brief Options for the expansion map.
162 : /// The outer boundary radius of the map is always set to
163 : /// the outer boundary of the Domain, so there is no option
164 : /// here to set the outer boundary radius.
165 1 : struct ExpansionMapOptions {
166 0 : using type = Options::Auto<ExpansionMapOptions, Options::AutoLabel::None>;
167 0 : static std::string name() { return "ExpansionMap"; }
168 0 : static constexpr Options::String help = {
169 : "Options for the expansion map. Specify 'None' to not use this map."};
170 0 : struct InitialValues {
171 0 : using type = std::array<double, 2>;
172 0 : static constexpr Options::String help = {
173 : "Initial value and deriv of expansion."};
174 : };
175 0 : struct AsymptoticVelocityOuterBoundary {
176 0 : using type = double;
177 0 : static constexpr Options::String help = {
178 : "The asymptotic velocity of the outer boundary."};
179 : };
180 0 : struct DecayTimescaleOuterBoundaryVelocity {
181 0 : using type = double;
182 0 : static constexpr Options::String help = {
183 : "The timescale for how fast the outer boundary velocity approaches "
184 : "its asymptotic value."};
185 : };
186 0 : using options = tmpl::list<InitialValues, AsymptoticVelocityOuterBoundary,
187 : DecayTimescaleOuterBoundaryVelocity>;
188 :
189 0 : std::array<double, 2> initial_values{
190 : std::numeric_limits<double>::signaling_NaN(),
191 : std::numeric_limits<double>::signaling_NaN()};
192 0 : double outer_boundary_velocity{
193 : std::numeric_limits<double>::signaling_NaN()};
194 0 : double outer_boundary_decay_time{
195 : std::numeric_limits<double>::signaling_NaN()};
196 : };
197 :
198 0 : struct RotationMapOptions {
199 0 : using type = Options::Auto<RotationMapOptions, Options::AutoLabel::None>;
200 0 : static std::string name() { return "RotationMap"; }
201 0 : static constexpr Options::String help = {
202 : "Options for a time-dependent rotation map about an arbitrary axis. "
203 : "Specify 'None' to not use this map."};
204 :
205 0 : struct InitialAngularVelocity {
206 0 : using type = std::array<double, 3>;
207 0 : static constexpr Options::String help = {"The initial angular velocity."};
208 : };
209 :
210 0 : using options = tmpl::list<InitialAngularVelocity>;
211 :
212 0 : std::array<double, 3> initial_angular_velocity{};
213 : };
214 :
215 : /// \brief Options for the Translation Map, the outer radius is always set to
216 : /// the outer boundary of the Domain, so there's no option needed for outer
217 : /// boundary.
218 1 : struct TranslationMapOptions {
219 0 : using type = Options::Auto<TranslationMapOptions, Options::AutoLabel::None>;
220 0 : static std::string name() { return "TranslationMap"; }
221 0 : static constexpr Options::String help = {
222 : "Options for a time-dependent translation map. Specify 'None' to not "
223 : "use this map."};
224 :
225 0 : struct InitialValues {
226 0 : using type = std::array<std::array<double, 3>, 3>;
227 0 : static constexpr Options::String help = {
228 : "Initial position, velocity and acceleration."};
229 : };
230 :
231 0 : using options = tmpl::list<InitialValues>;
232 :
233 0 : std::array<std::array<double, 3>, 3> initial_values{};
234 : };
235 :
236 : // We use a type alias here instead of defining the ShapeMapOptions struct
237 : // because there appears to be a bug in clang-10. If the definition of
238 : // ShapeMapOptions is here inside TimeDependentMapOptions, on clang-10 there
239 : // is a linking error that there is an undefined reference to
240 : // Options::Option::parse_as<TimeDependentMapOptions<A>> (and B). This doesn't
241 : // show up for GCC. If we put the definition of ShapeMapOptions outside of
242 : // TimeDependentMapOptions and just use a type alias here, the linking error
243 : // goes away.
244 : template <domain::ObjectLabel Object>
245 0 : using ShapeMapOptions =
246 : domain::creators::time_dependent_options::ShapeMapOptions<
247 : not IsCylindrical, Object>;
248 :
249 0 : using options =
250 : tmpl::list<InitialTime, ExpansionMapOptions, RotationMapOptions,
251 : TranslationMapOptions, ShapeMapOptions<domain::ObjectLabel::A>,
252 : ShapeMapOptions<domain::ObjectLabel::B>>;
253 0 : static constexpr Options::String help{
254 : "The options for all time dependent maps in a binary compact object "
255 : "domain. Specify 'None' to not use any time dependent maps."};
256 :
257 0 : TimeDependentMapOptions() = default;
258 :
259 0 : TimeDependentMapOptions(
260 : double initial_time,
261 : std::optional<ExpansionMapOptions> expansion_map_options,
262 : std::optional<RotationMapOptions> rotation_options,
263 : std::optional<TranslationMapOptions> translation_options,
264 : std::optional<ShapeMapOptions<domain::ObjectLabel::A>> shape_options_A,
265 : std::optional<ShapeMapOptions<domain::ObjectLabel::B>> shape_options_B,
266 : const Options::Context& context = {});
267 :
268 : /*!
269 : * \brief Create the function of time map using the options that were
270 : * provided to this class.
271 : *
272 : * Currently, this will add:
273 : *
274 : * - Expansion: `PiecewisePolynomial<2>`
275 : * - ExpansionOuterBoundary: `FixedSpeedCubic`
276 : * - Rotation: `QuaternionFunctionOfTime<3>`
277 : * - Translation: `PiecewisePolynomial<3>`
278 : * - SizeA/B: `PiecewisePolynomial<3>`
279 : * - ShapeA/B: `PiecewisePolynomial<2>`
280 : */
281 : std::unordered_map<std::string,
282 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
283 1 : create_functions_of_time(const std::unordered_map<std::string, double>&
284 : initial_expiration_times) const;
285 :
286 : /*!
287 : * \brief Construct the actual maps that will be used.
288 : *
289 : * Currently, this constructs a:
290 : *
291 : * - Rotation, Expansion, Translation: `RotScaleTrans<3>`
292 : * - ShapeA/B: `Shape` (with size FunctionOfTime)
293 : *
294 : * If the radii for an object are `std::nullopt`, this means that a Shape map
295 : * is not constructed for that object. An identity map will be used instead.
296 : * If \p IsCylindrical is true, only pass two radii for the inner/outer radius
297 : * of the object sphere. If it is false, pass three radii corresponding to the
298 : * excision radius, the outer radius of the inner sphere, and the radius of
299 : * the surrounding cube.
300 : */
301 1 : void build_maps(
302 : const std::array<std::array<double, 3>, 2>& centers,
303 : const std::optional<std::array<double, IsCylindrical ? 2 : 3>>&
304 : object_A_radii,
305 : const std::optional<std::array<double, IsCylindrical ? 2 : 3>>&
306 : object_B_radii,
307 : double envelope_radius, double domain_outer_radius);
308 :
309 : /*!
310 : * \brief Check whether options were specified in the constructor for the
311 : * shape map of this object
312 : */
313 1 : bool has_distorted_frame_options(domain::ObjectLabel object) const;
314 :
315 : /*!
316 : * \brief Type to pass to `grid_to_distorted_map()` and
317 : * `grid_to_inertial_map()`.
318 : *
319 : * \details If \p IsCylindrical is true, pass a `bool` for whether to include
320 : * a shape map or not. If it's false, pass a `std::optional<size_t>`. If this
321 : * has a value, then it includes a shape map. The `size_t` represents the
322 : * relative block number around each object in the BinaryCompactObject domain.
323 : * It should go from 0 to 11 for the 12 blocks surrounding each object.
324 : */
325 1 : using IncludeDistortedMapType =
326 : tmpl::conditional_t<IsCylindrical, bool, std::optional<size_t>>;
327 :
328 : /*!
329 : * \brief This will construct the map from `Frame::Distorted` to
330 : * `Frame::Inertial`
331 : *
332 : * If we are including a shape map, then this will be a `RotScaleTrans` map.
333 : * If we aren't, this returns a `nullptr`.
334 : *
335 : * \see IncludeDistortedMapType
336 : */
337 : template <domain::ObjectLabel Object>
338 1 : MapType<Frame::Distorted, Frame::Inertial> distorted_to_inertial_map(
339 : const IncludeDistortedMapType& include_distorted_map,
340 : bool use_rigid_map) const;
341 :
342 : /*!
343 : * \brief This will construct the maps from the `Frame::Grid` to the
344 : * `Frame::Distorted`.
345 : *
346 : * If we are including a shape map, then this will be a `Shape` map (with size
347 : * FunctionOfTime) for the templated `Object`. If we aren't, then this returns
348 : * a `nullptr`.
349 : *
350 : * \see IncludeDistortedMapType
351 : */
352 : template <domain::ObjectLabel Object>
353 1 : MapType<Frame::Grid, Frame::Distorted> grid_to_distorted_map(
354 : const IncludeDistortedMapType& include_distorted_map) const;
355 :
356 : /*!
357 : * \brief This will construct the entire map from the `Frame::Grid` to the
358 : * `Frame::Inertial`.
359 : *
360 : * If we are including a shape map, then this map will have a composition of a
361 : * `Shape` (with size FunctionOfTime) and `RotScaleTrans` map. If
362 : * not, there will only be a `RotScaleTrans` map.
363 : *
364 : * \see IncludeDistortedMapType
365 : */
366 : template <domain::ObjectLabel Object>
367 1 : MapType<Frame::Grid, Frame::Inertial> grid_to_inertial_map(
368 : const IncludeDistortedMapType& include_distorted_map,
369 : bool use_rigid_map) const;
370 :
371 : // Names are public because they need to be used when constructing maps in the
372 : // BCO domain creators themselves
373 0 : inline static const std::string expansion_name{"Expansion"};
374 0 : inline static const std::string expansion_outer_boundary_name{
375 : "ExpansionOuterBoundary"};
376 0 : inline static const std::string rotation_name{"Rotation"};
377 0 : inline static const std::string translation_name{"Translation"};
378 0 : inline static const std::array<std::string, 2> size_names{{"SizeA", "SizeB"}};
379 0 : inline static const std::array<std::string, 2> shape_names{
380 : {"ShapeA", "ShapeB"}};
381 :
382 : private:
383 0 : static size_t get_index(domain::ObjectLabel object);
384 :
385 0 : double initial_time_{std::numeric_limits<double>::signaling_NaN()};
386 0 : std::optional<ExpansionMapOptions> expansion_map_options_{};
387 0 : std::optional<RotationMapOptions> rotation_options_{};
388 0 : std::optional<TranslationMapOptions> translation_options_{};
389 0 : std::optional<ShapeMapOptions<domain::ObjectLabel::A>> shape_options_A_{};
390 0 : std::optional<ShapeMapOptions<domain::ObjectLabel::B>> shape_options_B_{};
391 0 : std::array<std::optional<double>, 2> inner_radii_{};
392 :
393 : // Maps
394 0 : std::optional<Expansion> expansion_map_{};
395 0 : std::optional<Rotation> rotation_map_{};
396 0 : std::optional<std::pair<RotScaleTrans, RotScaleTrans>> rot_scale_trans_map_{};
397 0 : using ShapeMapType =
398 : tmpl::conditional_t<IsCylindrical, std::array<std::optional<Shape>, 2>,
399 : std::array<std::array<std::optional<Shape>, 6>, 2>>;
400 0 : ShapeMapType shape_maps_{};
401 : };
402 : } // namespace domain::creators::bco
|
