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 <cstddef>
8 : #include <memory>
9 : #include <optional>
10 : #include <string>
11 : #include <unordered_map>
12 : #include <unordered_set>
13 : #include <variant>
14 : #include <vector>
15 :
16 : #include "DataStructures/Tensor/Tensor.hpp"
17 : #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
18 : #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp"
19 : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
20 : #include "Domain/CoordinateMaps/Identity.hpp"
21 : #include "Domain/Creators/DomainCreator.hpp"
22 : #include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp"
23 : #include "Domain/Domain.hpp"
24 : #include "Domain/Structure/DirectionMap.hpp"
25 : #include "Domain/Structure/ObjectLabel.hpp"
26 : #include "Options/Auto.hpp"
27 : #include "Options/Context.hpp"
28 : #include "Options/String.hpp"
29 : #include "Utilities/GetOutput.hpp"
30 : #include "Utilities/TMPL.hpp"
31 :
32 : /// \cond
33 : namespace domain {
34 : namespace CoordinateMaps {
35 : class Affine;
36 : enum class Distribution;
37 : template <size_t Dim>
38 : class Identity;
39 : class Interval;
40 : class PolarToCartesian;
41 : template <typename Map1, typename Map2>
42 : class ProductOf2Maps;
43 : template <typename Map1, typename Map2, typename Map3>
44 : class ProductOf3Maps;
45 : class SphericalToCartesianPfaffian;
46 : template <size_t VolumeDim>
47 : class Wedge;
48 : template <size_t VolumeDim>
49 : class DiscreteRotation;
50 : class UniformCylindricalEndcap;
51 : class UniformCylindricalFlatEndcap;
52 : class UniformCylindricalSide;
53 : } // namespace CoordinateMaps
54 :
55 : template <typename SourceFrame, typename TargetFrame, typename... Maps>
56 : class CoordinateMap;
57 :
58 : template <typename T>
59 : struct ExpandOverBlocks;
60 :
61 : namespace FunctionsOfTime {
62 : class FunctionOfTime;
63 : } // namespace FunctionsOfTime
64 : } // namespace domain
65 :
66 : namespace Frame {
67 : struct Grid;
68 : struct Distorted;
69 : struct Inertial;
70 : struct BlockLogical;
71 : } // namespace Frame
72 : /// \endcond
73 :
74 : namespace domain::creators {
75 :
76 : /*!
77 : * \ingroup ComputationalDomainGroup
78 : *
79 : * \brief A general domain for two compact objects based on cylinders.
80 : *
81 : * Creates a 3D Domain that represents a binary compact object
82 : * solution. This domain is described briefly in the Appendix of
83 : * \cite Buchman:2012dw, and is illustrated in Figure 20 of that
84 : * paper.
85 : *
86 : * In the code and options below, `ObjectA` and `ObjectB` refer to the
87 : * two compact objects. In the grid frame, `ObjectA` is located to the
88 : * right of (i.e. a more positive value of the x-coordinate than)
89 : * `ObjectB`. The inner edge of the Blocks surrounding each of
90 : * `ObjectA` and `ObjectB` is spherical in grid coordinates; the
91 : * user must specify the center and radius of this surface for both
92 : * `ObjectA` and `ObjectB`, and the user must specify the outer boundary
93 : * radius. The outer boundary is a sphere centered at the origin.
94 : *
95 : * This domain offers some grid anchors. See
96 : * `domain::creators::bco::create_grid_anchors` for which ones are offered.
97 : *
98 : * Note that Figure 20 of \cite Buchman:2012dw illustrates additional
99 : * spherical shells inside the "EA" and "EB" blocks, and the caption
100 : * of Figure 20 indicates that there are additional spherical shells
101 : * outside the "CA" and "CB" blocks; `CylindricalBinaryCompactObject`
102 : * has these extra shells inside "EA" only if the option `InnerSphereA`
103 : * is not `None`, and it has the extra shells inside "EB" only if the option
104 : * `InnerSphereB` is not `None`. If the shells are absent, then the "EA" and
105 : * "EB" blocks extend to the excision boundaries.
106 : *
107 : * The Blocks are named as follows:
108 : * - Each of CAFilledCylinder, EAFilledCylinder, EBFilledCylinder,
109 : * MAFilledCylinder, MBFilledCylinder, and CBFilledCylinder are filled
110 : * cylindrical endcaps made of a single cylindrical block.
111 : * - Each of CACylinder, EACylinder, EBCylinder, and CBCylinder are hollow
112 : * ylindrical shells made of a single cylindrical block.
113 : * - The Block group called "Outer" consists of all the CA and CB blocks.
114 : * - OuterShell0 is the single shell in a Block group called "OuterSphere" and
115 : * it borders the outer boundary.
116 : * - The Block group called "InnerA" consists of all the EA, and MA
117 : * blocks. They all border the inner boundary "A" if `InnerSphereA` is not
118 : * `None`.
119 : * - If `InnerSphereA` is not `None`, InnerAShell0 is the single shell in a
120 : * Block group called "InnerSphereA" and it borders the inner excision
121 : * boundary "A".
122 : * - The Block group called "InnerB" consists of all the EB, and MB
123 : * blocks. They all border the inner boundary "B" if `InnerSphereB` is not
124 : * `None`.
125 : * - If `InnerSphereB` is true, InnerBShell0 is the single shell in a
126 : * Block group called "InnerSphereB" and it borders the inner excision
127 : * boundary "B".
128 : *
129 : * If \f$c_A\f$ and \f$c_B\f$ are the input parameters center_A and
130 : * center_B, \f$r_A\f$ and \f$r_B\f$ are the input parameters radius_A and
131 : * radius_B, and \f$R\f$ is the outer boundary radius, we demand the
132 : * following restrictions on parameters:
133 : * - \f$c_A^0>0\f$; this is a convention to simplify the code.
134 : * - \f$c_B^0<0\f$; this is a convention to simplify the code.
135 : * - \f$|c_A^0|\le|c_B^0|\f$. We should roughly have \f$r_A c_A^0 + r_B c_B^0\f$
136 : * close to zero; that is, for BBHs (where \f$r_A\f$ is roughly twice the
137 : * mass of the heavier object A, and \f$r_B\f$ is roughly twice the mass
138 : * of the lighter object B) the center of mass should be roughly
139 : * at the origin.
140 : * - \f$0 < r_B < r_A\f$
141 : * - \f$R \ge 3(|c_A^0|-|c_B^0|)\f$; otherwise the blocks will be too compressed
142 : * near the outer boundary.
143 : *
144 : * All time dependent maps are optional to specify. To include a map, specify
145 : * its options. Otherwise specify `None` for that map. You can also turn off
146 : * time dependent maps all together by specifying `None` for the
147 : * `TimeDependentMaps` option. See
148 : * `domain::creators::bco::TimeDependentMapOptions`. This class must pass a
149 : * template parameter of `true` to
150 : * `domain::creators::bco::TimeDependentMapOptions`.
151 : */
152 1 : class CylindricalBinaryCompactObject : public DomainCreator<3> {
153 : public:
154 0 : using unit_cylinder_map =
155 : CoordinateMaps::ProductOf3Maps<CoordinateMaps::Affine,
156 : CoordinateMaps::Identity<1>,
157 : CoordinateMaps::Interval>;
158 0 : using polar_to_cartesian_map =
159 : CoordinateMaps::ProductOf2Maps<CoordinateMaps::PolarToCartesian,
160 : CoordinateMaps::Identity<1>>;
161 :
162 0 : using maps_list = tmpl::flatten<tmpl::list<
163 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
164 : unit_cylinder_map, polar_to_cartesian_map,
165 : CoordinateMaps::DiscreteRotation<3>,
166 : CoordinateMaps::UniformCylindricalEndcap,
167 : CoordinateMaps::DiscreteRotation<3>>,
168 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
169 : unit_cylinder_map, polar_to_cartesian_map,
170 : CoordinateMaps::DiscreteRotation<3>,
171 : CoordinateMaps::UniformCylindricalFlatEndcap,
172 : CoordinateMaps::DiscreteRotation<3>>,
173 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
174 : unit_cylinder_map, polar_to_cartesian_map,
175 : CoordinateMaps::DiscreteRotation<3>,
176 : CoordinateMaps::UniformCylindricalSide,
177 : CoordinateMaps::DiscreteRotation<3>>,
178 : domain::CoordinateMap<
179 : Frame::BlockLogical, Frame::Inertial,
180 : domain::CoordinateMaps::ProductOf2Maps<CoordinateMaps::Interval,
181 : CoordinateMaps::Identity<2>>,
182 : domain::CoordinateMaps::SphericalToCartesianPfaffian,
183 : CoordinateMaps::ProductOf3Maps<CoordinateMaps::Affine,
184 : CoordinateMaps::Affine,
185 : CoordinateMaps::Affine>>,
186 : bco::TimeDependentMapOptions<true>::maps_list>>;
187 :
188 0 : struct CenterA {
189 0 : using type = std::array<double, 3>;
190 0 : static constexpr Options::String help = {
191 : "Grid coordinates of center for Object A, which is at x>0."};
192 : };
193 0 : struct CenterB {
194 0 : using type = std::array<double, 3>;
195 0 : static constexpr Options::String help = {
196 : "Grid coordinates of center for Object B, which is at x<0."};
197 : };
198 0 : struct RadiusA {
199 0 : using type = double;
200 0 : static constexpr Options::String help = {
201 : "Grid-coordinate radius of grid boundary around Object A."};
202 : };
203 0 : struct RadiusB {
204 0 : using type = double;
205 0 : static constexpr Options::String help = {
206 : "Grid-coordinate radius of grid boundary around Object B."};
207 : };
208 :
209 : /// Options for InnerSphereA and InnerSphereB
210 1 : struct InnerSphere {
211 0 : static constexpr Options::String help = {
212 : "Options for an inner sphere, an extra spherical layer of Blocks "
213 : "around one of the objects. Specify 'None' to exclude."};
214 :
215 0 : struct OuterRadius {
216 0 : using type = Options::Auto<double>;
217 0 : static constexpr Options::String help = {
218 : "Outer radius of the inner sphere, or Auto to have it chosen for "
219 : "you."};
220 : };
221 :
222 0 : struct RadialDistribution {
223 0 : using type = domain::CoordinateMaps::Distribution;
224 0 : static constexpr Options::String help = {
225 : "The radial distribution of grid points for the inner sphere."};
226 : };
227 :
228 0 : using options = tmpl::list<OuterRadius, RadialDistribution>;
229 :
230 0 : InnerSphere() = default;
231 0 : explicit InnerSphere(
232 : std::optional<double> outer_radius,
233 : const domain::CoordinateMaps::Distribution radial_distribution)
234 : : outer_radius_(outer_radius),
235 : radial_distribution_(radial_distribution) {}
236 :
237 0 : std::optional<double> outer_radius_;
238 0 : domain::CoordinateMaps::Distribution radial_distribution_{
239 : domain::CoordinateMaps::Distribution::Logarithmic};
240 : };
241 :
242 0 : struct InnerSphereA {
243 0 : using type = Options::Auto<InnerSphere, Options::AutoLabel::None>;
244 0 : static constexpr Options::String help = {"Options for InnerSphereA."};
245 : };
246 :
247 0 : struct InnerSphereB {
248 0 : using type = Options::Auto<InnerSphere, Options::AutoLabel::None>;
249 0 : static constexpr Options::String help = {"Options for InnerSphereB."};
250 : };
251 :
252 0 : struct OuterSphereOptions {
253 : public:
254 0 : static constexpr Options::String help = {"Options for OuterSphere."};
255 :
256 0 : struct OuterRadius {
257 0 : using type = double;
258 0 : static constexpr Options::String help = {"Outer radius of OuterSphere."};
259 : };
260 :
261 0 : using options = tmpl::list<OuterRadius>;
262 :
263 0 : OuterSphereOptions() = default;
264 0 : explicit OuterSphereOptions(const double outer_radius)
265 : : outer_radius_(outer_radius) {}
266 :
267 0 : double outer_radius_;
268 : };
269 :
270 0 : struct OuterSphere {
271 0 : using type = OuterSphereOptions;
272 0 : static constexpr Options::String help = {"Options for OuterSphere."};
273 : };
274 :
275 0 : struct InitialRefinement {
276 0 : using type = std::variant<
277 : size_t, std::unordered_map<
278 : std::string, std::variant<std::array<size_t, 2>, size_t>>>;
279 0 : static constexpr Options::String help = {
280 : "Initial refinement level. Specify one of: a single number for global "
281 : "refinement or a list of entries for every block/group in the domain, "
282 : "where spherical shell block/group entries are specified with one "
283 : "number and cylindrical block/group entries are speicfied with an "
284 : "array<2> representing [radial refinement, z refinement]. \n\nNote "
285 : "that the z direction in cylinder blocks will roughly correspond to "
286 : "refinement in a direction parallel to the axis of separation between "
287 : "the two objects. Because filled cylinder blocks lie along the axis of "
288 : "separation but hollow cylinder blocks wrap around it, refinement in z "
289 : "leads to refinement in different spherical coordinate directions in "
290 : "the global spherical coordinates. More specifically, z refinement in "
291 : "filled cylinders (e.g. EAFilledCylinder) will roughly correspond to "
292 : "radial refinement in global spherical coordinates, but in hollow "
293 : "cylinders, it will behave more like angular refinement in global "
294 : "spherical coordinates that is perpendicular to the cylinder's local "
295 : "angular direction."};
296 : };
297 :
298 0 : struct InitialGridPoints {
299 0 : using type = std::variant<
300 : size_t,
301 : std::unordered_map<std::string, std::variant<std::array<size_t, 3>,
302 : std::array<size_t, 2>>>>;
303 0 : static constexpr Options::String help = {
304 : "Initial number of grid points. Specify one of the following:"
305 : "\n\t- a single number"
306 : "\n\t- lists for blocks and/or block groups as follows:"
307 : "\n\t\t- [r, l_max] for spherical shell blocks and groups"
308 : "\n\t\t- [r, theta, z] for hollow cylinder blocks, where theta must be "
309 : "\n\t\t odd"
310 : "\n\t\t- [r, z] or [r, theta, z] for filled cylinder blocks and groups "
311 : "\n\t\t containing them, where r must be > 2, theta must be odd, and "
312 : "\n\t\t (4 * r - 3) >= theta\n\n"
313 : "While the most verbose, the best choice for a production run is "
314 : "likely to specify a list for each cylindrical block instead of each "
315 : "cylindrical block group. If you set a whole group (e.g. InnerA) using "
316 : "[r, z], the interfaces between its filled cylinders "
317 : "(e.g. EAFilledCylinder) and its hollow cylinders (e.g. EACylinder) "
318 : "may not have similar resolution on each side unless r and z are "
319 : "close. This is because at these interfaces, the z direction in filled "
320 : "cylinders lines up with the radial direction in hollow cylinders. To "
321 : "get the resolution on either side of these interfaces to match well, "
322 : "you either want to set a cylindrical block group with r and z close "
323 : "in value or set the individual cylindrical blocks for more freedom. "
324 : "Also note that any h refinement in these cylindrical blocks will also "
325 : "be in similarly different directions at the interface, which affects "
326 : "this picture of trying to match the p refinement at the interface of "
327 : "hollow and filled cylinders."};
328 : };
329 :
330 0 : struct BoundaryConditions {
331 0 : static constexpr Options::String help = "The boundary conditions to apply.";
332 : };
333 : template <typename BoundaryConditionsBase>
334 0 : struct InnerBoundaryCondition {
335 0 : static std::string name() { return "InnerBoundary"; }
336 0 : static constexpr Options::String help =
337 : "Options for the inner boundary conditions.";
338 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
339 0 : using group = BoundaryConditions;
340 : };
341 :
342 : template <typename BoundaryConditionsBase>
343 0 : struct OuterBoundaryCondition {
344 0 : static std::string name() { return "OuterBoundary"; }
345 0 : static constexpr Options::String help =
346 : "Options for the outer boundary conditions.";
347 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
348 0 : using group = BoundaryConditions;
349 : };
350 :
351 0 : struct TimeDependentMaps {
352 0 : using type = Options::Auto<bco::TimeDependentMapOptions<true>,
353 : Options::AutoLabel::None>;
354 0 : static constexpr Options::String help =
355 : bco::TimeDependentMapOptions<true>::help;
356 : };
357 :
358 : template <typename Metavariables>
359 0 : using options = tmpl::append<
360 : tmpl::list<CenterA, CenterB, RadiusA, RadiusB, InitialRefinement,
361 : InitialGridPoints, OuterSphere, InnerSphereA, InnerSphereB,
362 : TimeDependentMaps>,
363 : tmpl::conditional_t<
364 : domain::BoundaryConditions::has_boundary_conditions_base_v<
365 : typename Metavariables::system>,
366 : tmpl::list<
367 : InnerBoundaryCondition<
368 : domain::BoundaryConditions::get_boundary_conditions_base<
369 : typename Metavariables::system>>,
370 : OuterBoundaryCondition<
371 : domain::BoundaryConditions::get_boundary_conditions_base<
372 : typename Metavariables::system>>>,
373 : tmpl::list<>>>;
374 :
375 0 : static constexpr Options::String help{
376 : "The CylindricalBinaryCompactObject domain is a general domain for "
377 : "two compact objects. The user must provide the (grid-frame) "
378 : "centers and radii of the spherical inner edge of the grid surrounding "
379 : "each of the two compact objects A and B."};
380 :
381 0 : CylindricalBinaryCompactObject(
382 : std::array<double, 3> center_A, std::array<double, 3> center_B,
383 : double radius_A, double radius_B,
384 : const typename InitialRefinement::type& initial_refinement,
385 : const typename InitialGridPoints::type& initial_grid_points,
386 : typename OuterSphere::type outer_sphere,
387 : std::optional<InnerSphere> inner_sphere_A = std::nullopt,
388 : std::optional<InnerSphere> inner_sphere_B = std::nullopt,
389 : std::optional<bco::TimeDependentMapOptions<true>> time_dependent_options =
390 : std::nullopt,
391 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
392 : inner_boundary_condition = nullptr,
393 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
394 : outer_boundary_condition = nullptr,
395 : const Options::Context& context = {});
396 :
397 0 : CylindricalBinaryCompactObject() = default;
398 0 : CylindricalBinaryCompactObject(const CylindricalBinaryCompactObject&) =
399 : delete;
400 0 : CylindricalBinaryCompactObject(CylindricalBinaryCompactObject&&) = default;
401 0 : CylindricalBinaryCompactObject& operator=(
402 : const CylindricalBinaryCompactObject&) = delete;
403 0 : CylindricalBinaryCompactObject& operator=(CylindricalBinaryCompactObject&&) =
404 : default;
405 0 : ~CylindricalBinaryCompactObject() override = default;
406 :
407 0 : Domain<3> create_domain() const override;
408 :
409 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
410 1 : grid_anchors() const override {
411 : return grid_anchors_;
412 : }
413 :
414 : std::vector<DirectionMap<
415 : 3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
416 1 : external_boundary_conditions() const override;
417 :
418 1 : std::vector<std::array<size_t, 3>> initial_extents() const override;
419 :
420 1 : std::vector<std::array<size_t, 3>> initial_refinement_levels() const override;
421 :
422 1 : auto functions_of_time(const std::unordered_map<std::string, double>&
423 : initial_expiration_times = {}) const
424 : -> std::unordered_map<
425 : std::string,
426 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>> override;
427 :
428 1 : std::vector<std::string> block_names() const override { return block_names_; }
429 :
430 : std::unordered_map<std::string, std::unordered_set<std::string>>
431 1 : block_groups() const override {
432 : return block_groups_;
433 : }
434 :
435 : private:
436 : // Note that center_A_ and center_B_ are rotated with respect to the
437 : // input centers (which are in the grid frame), so that we can
438 : // construct the map in a frame where the centers are offset in the
439 : // z direction. At the end, there will be another rotation back to
440 : // the grid frame (where the centers are offset in the x direction).
441 0 : std::array<double, 3> center_A_{};
442 0 : std::array<double, 3> center_B_{};
443 0 : double radius_A_{};
444 0 : double radius_B_{};
445 0 : double outer_radius_A_{};
446 0 : double outer_radius_B_{};
447 0 : typename std::vector<std::array<size_t, 3>> initial_refinement_{};
448 0 : typename std::vector<std::array<size_t, 3>> initial_grid_points_{};
449 0 : typename OuterSphere::type outer_sphere_{};
450 0 : std::optional<InnerSphere> inner_sphere_A_{};
451 0 : std::optional<InnerSphere> inner_sphere_B_{};
452 : // cut_spheres_offset_factor_ is eta in Eq. (A.9) of
453 : // https://arxiv.org/abs/1206.3015. cut_spheres_offset_factor_
454 : // could be set to unity to simplify the equations. Here we fix it
455 : // to the value 0.99 used in SpEC, so that we reproduce SpEC's
456 : // domain decomposition.
457 0 : double cut_spheres_offset_factor_{0.99};
458 : // z_cutting_plane_ is x_C in Eq. (A.9) of
459 : // https://arxiv.org/abs/1206.3015 (but rotated to the z-axis).
460 0 : double z_cutting_plane_{};
461 0 : size_t number_of_blocks_{};
462 0 : std::unordered_map<std::string, size_t> block_positions_{};
463 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
464 0 : inner_boundary_condition_;
465 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
466 0 : outer_boundary_condition_;
467 0 : std::vector<std::string> block_names_{};
468 : std::unordered_map<std::string, std::unordered_set<std::string>>
469 0 : block_groups_{};
470 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
471 0 : grid_anchors_{};
472 : // FunctionsOfTime options
473 0 : std::optional<bco::TimeDependentMapOptions<true>> time_dependent_options_{};
474 : };
475 : } // namespace domain::creators
|
