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