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/TimeDependentOptions/ShapeMap.hpp"
21 : #include "Domain/Creators/TimeDependentOptions/Sphere.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 : namespace domain::creators::bco {
42 : namespace detail {
43 : // Convenience type alias
44 : template <typename... Maps>
45 : using gi_map = domain::CoordinateMap<Frame::Grid, Frame::Inertial, Maps...>;
46 : template <typename... Maps>
47 : using gd_map = domain::CoordinateMap<Frame::Grid, Frame::Distorted, Maps...>;
48 : template <typename... Maps>
49 : using di_map =
50 : domain::CoordinateMap<Frame::Distorted, Frame::Inertial, Maps...>;
51 :
52 : template <typename List>
53 : struct power_set {
54 : using rest = typename power_set<tmpl::pop_front<List>>::type;
55 : using type = tmpl::append<
56 : rest, tmpl::transform<rest, tmpl::lazy::push_front<
57 : tmpl::_1, tmpl::pin<tmpl::front<List>>>>>;
58 : };
59 :
60 : template <>
61 : struct power_set<tmpl::list<>> {
62 : using type = tmpl::list<tmpl::list<>>;
63 : };
64 :
65 : template <typename SourceFrame, typename TargetFrame, typename Maps>
66 : using produce_all_maps_helper =
67 : tmpl::wrap<tmpl::push_front<Maps, SourceFrame, TargetFrame>,
68 : domain::CoordinateMap>;
69 :
70 : /*
71 : * This will produce all the possible combinations of maps we will need for the
72 : * maps_list. It does so using a power set. Say you have 3 maps, Map1, Map2,
73 : * Map3. The end result of this will be a list with the following map
74 : * combinations where these combinations may not appear in the same position in
75 : * the list, but the order of each combination is fixed:
76 : * - Map1
77 : * - Map2
78 : * - Map3
79 : * - Map1, Map2
80 : * - Map1, Map3
81 : * - Map2, Map3
82 : * - Map1, Map2, Map3
83 : */
84 : template <typename SourceFrame, typename TargetFrame, typename... Maps>
85 : using produce_all_maps = tmpl::transform<
86 : tmpl::remove<typename power_set<tmpl::list<Maps...>>::type, tmpl::list<>>,
87 : tmpl::bind<produce_all_maps_helper, tmpl::pin<SourceFrame>,
88 : tmpl::pin<TargetFrame>, tmpl::_1>>;
89 : } // namespace detail
90 :
91 : /*!
92 : * \brief This holds all options related to the time dependent maps of the
93 : * binary compact object domains.
94 : *
95 : * \details Since both domains will have the same (overall) time dependent maps,
96 : * their options are going to be the same as well. The options won't be exactly
97 : * the same though, so there is a \p IsCylindrical template parameter to
98 : * distinguish.
99 : *
100 : * This class will create the FunctionsOfTime needed for the binary compact
101 : * object domains as well as the actual CoordinateMap%s themselves
102 : *
103 : * \note This struct contains no information about what blocks the time
104 : * dependent maps will go in.
105 : */
106 : template <bool IsCylindrical>
107 1 : struct TimeDependentMapOptions {
108 : private:
109 : template <typename SourceFrame, typename TargetFrame>
110 0 : using MapType =
111 : std::unique_ptr<domain::CoordinateMapBase<SourceFrame, TargetFrame, 3>>;
112 : // Time-dependent maps
113 0 : using Expansion = domain::CoordinateMaps::TimeDependent::CubicScale<3>;
114 0 : using Rotation = domain::CoordinateMaps::TimeDependent::Rotation<3>;
115 0 : using RotScaleTrans = domain::CoordinateMaps::TimeDependent::RotScaleTrans<3>;
116 0 : using Shape = domain::CoordinateMaps::TimeDependent::Shape;
117 0 : using Identity = domain::CoordinateMaps::Identity<3>;
118 :
119 : public:
120 0 : using maps_list = tmpl::append<
121 : // We need this odd-one-out Identity map because all maps are optional to
122 : // specify. It's possible to specify a shape map, but no expansion or
123 : // rotation map. So we have a grid to distorted map, and need an identity
124 : // distorted to inertial map. We don't need an identity grid to distorted
125 : // map because if a user requests the grid to distorted frame map with a
126 : // distorted frame, but didn't specify shape map options, an error occurs.
127 : tmpl::list<detail::di_map<Identity>>,
128 : detail::produce_all_maps<Frame::Grid, Frame::Inertial, Shape, Expansion,
129 : Rotation>,
130 : detail::produce_all_maps<Frame::Grid, Frame::Distorted, Shape>,
131 : detail::produce_all_maps<Frame::Distorted, Frame::Inertial, Expansion,
132 : Rotation>,
133 : detail::produce_all_maps<Frame::Grid, Frame::Inertial, Shape,
134 : RotScaleTrans>,
135 : detail::produce_all_maps<Frame::Distorted, Frame::Inertial,
136 : RotScaleTrans>>;
137 :
138 : /// \brief The initial time of the functions of time.
139 1 : struct InitialTime {
140 0 : using type = double;
141 0 : static constexpr Options::String help = {
142 : "The initial time of the functions of time"};
143 : };
144 :
145 : /// \brief Options for the expansion map.
146 : /// The outer boundary radius of the map is always set to
147 : /// the outer boundary of the Domain, so there is no option
148 : /// here to set the outer boundary radius.
149 1 : struct ExpansionMapOptions {
150 0 : using type = Options::Auto<ExpansionMapOptions, Options::AutoLabel::None>;
151 0 : static std::string name() { return "ExpansionMap"; }
152 0 : static constexpr Options::String help = {
153 : "Options for the expansion map. Specify 'None' to not use this map."};
154 0 : struct InitialValues {
155 0 : using type = std::array<double, 2>;
156 0 : static constexpr Options::String help = {
157 : "Initial value and deriv of expansion."};
158 : };
159 0 : struct AsymptoticVelocityOuterBoundary {
160 0 : using type = double;
161 0 : static constexpr Options::String help = {
162 : "The asymptotic velocity of the outer boundary."};
163 : };
164 0 : struct DecayTimescaleOuterBoundaryVelocity {
165 0 : using type = double;
166 0 : static constexpr Options::String help = {
167 : "The timescale for how fast the outer boundary velocity approaches "
168 : "its asymptotic value."};
169 : };
170 0 : using options = tmpl::list<InitialValues, AsymptoticVelocityOuterBoundary,
171 : DecayTimescaleOuterBoundaryVelocity>;
172 0 : ExpansionMapOptions() = default;
173 0 : ExpansionMapOptions(std::array<double, 2> initial_values_in,
174 : double outer_boundary_velocity_in,
175 : double outer_boundary_decay_time_in)
176 : : initial_values(initial_values_in),
177 : outer_boundary_velocity(outer_boundary_velocity_in),
178 : outer_boundary_decay_time(outer_boundary_decay_time_in) {}
179 :
180 0 : std::array<double, 2> initial_values{
181 : std::numeric_limits<double>::signaling_NaN(),
182 : std::numeric_limits<double>::signaling_NaN()};
183 0 : double outer_boundary_velocity{
184 : std::numeric_limits<double>::signaling_NaN()};
185 0 : double outer_boundary_decay_time{
186 : std::numeric_limits<double>::signaling_NaN()};
187 : };
188 :
189 0 : struct RotationMapOptions {
190 0 : using type = Options::Auto<RotationMapOptions, Options::AutoLabel::None>;
191 0 : static std::string name() { return "RotationMap"; }
192 0 : static constexpr Options::String help = {
193 : "Options for a time-dependent rotation map about an arbitrary axis. "
194 : "Specify 'None' to not use this map."};
195 :
196 0 : struct InitialAngularVelocity {
197 0 : using type = std::array<double, 3>;
198 0 : static constexpr Options::String help = {"The initial angular velocity."};
199 : };
200 :
201 0 : using options = tmpl::list<InitialAngularVelocity>;
202 :
203 0 : RotationMapOptions() = default;
204 0 : explicit RotationMapOptions(
205 : std::array<double, 3> initial_angular_velocity_in)
206 : : initial_angular_velocity(initial_angular_velocity_in) {}
207 :
208 0 : std::array<double, 3> initial_angular_velocity{};
209 : };
210 :
211 : /// \brief Options for the Translation Map, the outer radius is always set to
212 : /// the outer boundary of the Domain, so there's no option needed for outer
213 : /// boundary.
214 1 : struct TranslationMapOptions {
215 0 : using type = Options::Auto<TranslationMapOptions, Options::AutoLabel::None>;
216 0 : static std::string name() { return "TranslationMap"; }
217 0 : static constexpr Options::String help = {
218 : "Options for a time-dependent translation map. Specify 'None' to not "
219 : "use this map."};
220 :
221 0 : struct InitialValues {
222 0 : using type = std::array<std::array<double, 3>, 3>;
223 0 : static constexpr Options::String help = {
224 : "Initial position, velocity and acceleration."};
225 : };
226 :
227 0 : using options = tmpl::list<InitialValues>;
228 0 : TranslationMapOptions() = default;
229 0 : explicit TranslationMapOptions(
230 : std::array<std::array<double, 3>, 3> initial_values_in)
231 : : initial_values(initial_values_in) {}
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_map_options,
263 : std::optional<TranslationMapOptions> translation_map_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 : * When `UseWorldtube` is set to true, they are
282 : *
283 : * - Expansion: `IntegratedFunctionOfTime`
284 : * - ExpansionOuterBoundary: `FixedSpeedCubic`
285 : * - Rotation: `IntegratedFunctionOfTime`
286 : * - Translation: None
287 : * - SizeA/B: IntegratedFunctionOfTime
288 : * - ShapeA/B: `PiecewisePolynomial<2>`
289 : */
290 : template <bool UseWorldtube = false>
291 : std::unordered_map<std::string,
292 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
293 1 : create_functions_of_time(const std::unordered_map<std::string, double>&
294 : initial_expiration_times) const;
295 :
296 : /*!
297 : * \brief Construct the actual maps that will be used.
298 : *
299 : * Currently, this constructs a:
300 : *
301 : * - Rotation, Expansion, Translation: `RotScaleTrans<3>`
302 : * - ShapeA/B: `Shape` (with size FunctionOfTime)
303 : *
304 : * If the radii for an object are `std::nullopt`, this means that a Shape map
305 : * is not constructed for that object. An identity map will be used instead.
306 : * This happens when the object is covered by a Cartesian cube, rather than a
307 : * sphere.
308 : * If \p IsCylindrical is true, only pass two radii for the inner/outer radius
309 : * of the object sphere. If it is false, pass three radii corresponding to the
310 : * excision radius, the outer radius of the inner sphere, and the radius of
311 : * the surrounding cube.
312 : *
313 : * If a shape map is requested and the object is excised, then the shape map
314 : * will deform the inner excision surface of the object. This deformation
315 : * either extends to the outer radius of the sphere, or further to the edge of
316 : * the cube that surrounds the sphere, depending on the `TransitionEndsAtCube`
317 : * option.
318 : * If the object is filled, then the shape map will deform the outer radius of
319 : * the sphere that covers the object and the `TransitionEndsAtCube` option
320 : * must be `true`.
321 : */
322 1 : void build_maps(
323 : const std::array<std::array<double, 3>, 2>& object_centers,
324 : const std::optional<std::array<double, 3>>& cube_A_center,
325 : const std::optional<std::array<double, 3>>& cube_B_center,
326 : const std::optional<std::array<double, IsCylindrical ? 2 : 3>>&
327 : object_A_radii,
328 : const std::optional<std::array<double, IsCylindrical ? 2 : 3>>&
329 : object_B_radii,
330 : bool object_A_filled, bool object_B_filled, double envelope_radius,
331 : double domain_outer_radius);
332 :
333 : /*!
334 : * \brief Check whether options were specified in the constructor for the
335 : * shape map of this object
336 : */
337 1 : bool has_distorted_frame_options(domain::ObjectLabel object) const;
338 :
339 : /*!
340 : * \brief Type to pass to `grid_to_distorted_map()` and
341 : * `grid_to_inertial_map()`.
342 : *
343 : * \details If \p IsCylindrical is true, pass a `bool` for whether to include
344 : * a shape map or not. If it's false, pass a `std::optional<size_t>`. If this
345 : * has a value, then it includes a shape map. The `size_t` represents the
346 : * relative block number around each object in the BinaryCompactObject domain.
347 : * It should go from 0 to 11 for the 12 blocks surrounding each object.
348 : */
349 1 : using IncludeDistortedMapType =
350 : tmpl::conditional_t<IsCylindrical, bool, std::optional<size_t>>;
351 :
352 : /*!
353 : * \brief This will construct the map from `Frame::Distorted` to
354 : * `Frame::Inertial`
355 : *
356 : * If we are including a shape map, then this will be a `RotScaleTrans` map.
357 : * If we aren't, this returns a `nullptr`.
358 : *
359 : * \see IncludeDistortedMapType
360 : */
361 : template <domain::ObjectLabel Object>
362 1 : MapType<Frame::Distorted, Frame::Inertial> distorted_to_inertial_map(
363 : const IncludeDistortedMapType& include_distorted_map,
364 : bool use_rigid_map) const;
365 :
366 : /*!
367 : * \brief This will construct the maps from the `Frame::Grid` to the
368 : * `Frame::Distorted`.
369 : *
370 : * If we are including a shape map, then this will be a `Shape` map (with size
371 : * FunctionOfTime) for the templated `Object`. If we aren't, then this returns
372 : * a `nullptr`.
373 : *
374 : * \see IncludeDistortedMapType
375 : */
376 : template <domain::ObjectLabel Object>
377 1 : MapType<Frame::Grid, Frame::Distorted> grid_to_distorted_map(
378 : const IncludeDistortedMapType& include_distorted_map) const;
379 :
380 : /*!
381 : * \brief This will construct the entire map from the `Frame::Grid` to the
382 : * `Frame::Inertial`.
383 : *
384 : * If we are including a shape map, then this map will have a composition of a
385 : * `Shape` (with size FunctionOfTime) and `RotScaleTrans` map. If
386 : * not, there will only be a `RotScaleTrans` map.
387 : *
388 : * \see IncludeDistortedMapType
389 : */
390 : template <domain::ObjectLabel Object>
391 1 : MapType<Frame::Grid, Frame::Inertial> grid_to_inertial_map(
392 : const IncludeDistortedMapType& include_distorted_map,
393 : bool use_rigid_map) const;
394 :
395 : // Names are public because they need to be used when constructing maps in the
396 : // BCO domain creators themselves
397 0 : inline static const std::string expansion_name{"Expansion"};
398 0 : inline static const std::string expansion_outer_boundary_name{
399 : "ExpansionOuterBoundary"};
400 0 : inline static const std::string rotation_name{"Rotation"};
401 0 : inline static const std::string translation_name{"Translation"};
402 0 : inline static const std::array<std::string, 2> size_names{{"SizeA", "SizeB"}};
403 0 : inline static const std::array<std::string, 2> shape_names{
404 : {"ShapeA", "ShapeB"}};
405 :
406 : private:
407 0 : static size_t get_index(domain::ObjectLabel object);
408 :
409 0 : double initial_time_{std::numeric_limits<double>::signaling_NaN()};
410 0 : std::optional<ExpansionMapOptions> expansion_map_options_{};
411 0 : std::optional<RotationMapOptions> rotation_map_options_{};
412 0 : std::optional<TranslationMapOptions> translation_map_options_{};
413 0 : std::optional<ShapeMapOptions<domain::ObjectLabel::A>> shape_options_A_{};
414 0 : std::optional<ShapeMapOptions<domain::ObjectLabel::B>> shape_options_B_{};
415 0 : std::array<std::optional<double>, 2> inner_radii_{};
416 :
417 : // Maps
418 0 : std::optional<Expansion> expansion_map_{};
419 0 : std::optional<Rotation> rotation_map_{};
420 0 : std::optional<std::pair<RotScaleTrans, RotScaleTrans>> rot_scale_trans_map_{};
421 0 : using ShapeMapType =
422 : tmpl::conditional_t<IsCylindrical, std::array<std::optional<Shape>, 2>,
423 : std::array<std::array<std::optional<Shape>, 12>, 2>>;
424 0 : ShapeMapType shape_maps_{};
425 :
426 : // helper function that creates the functions of time used by the worldtube
427 : std::unordered_map<std::string,
428 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
429 0 : create_worldtube_functions_of_time() const;
430 : };
431 : } // namespace domain::creators::bco
|
