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<size_t, std::unordered_map<std::string, size_t>>;
277 0 : static constexpr Options::String help = {
278 : "Initial refinement level. Specify one of: a single number or a list "
279 : "of single numbers for every block group in the domain, every block "
280 : "name in the domain, or a mix of block groups and blocks. Each single "
281 : "number represents the radial refinement for spherical shell blocks "
282 : "and z refinement for cylindrical blocks.\n\nNote that the z direction "
283 : "in cylinder blocks will roughly correspond to refinement in a "
284 : "direction parallel to the axis of separation between the two objects. "
285 : "Because filled cylinder blocks lie along the axis of separation but "
286 : "hollow cylinder blocks wrap around it, refinement in z leads to "
287 : "refinement in different spherical coordinate directions in the "
288 : "global spherical coordinates. More specifically, z refinement in "
289 : "filled cylinders (e.g. EAFilledCylinder) will roughly correspond to "
290 : "radial refinement in global spherical coordinates, but in hollow "
291 : "cylinders, it will behave more like angular refinement in global "
292 : "spherical coordinates that is perpendicular to the cylinder's local "
293 : "angular direction."};
294 : };
295 :
296 0 : struct InitialGridPoints {
297 0 : using type = std::variant<
298 : size_t,
299 : std::unordered_map<std::string, std::variant<std::array<size_t, 3>,
300 : std::array<size_t, 2>>>>;
301 0 : static constexpr Options::String help = {
302 : "Initial number of grid points. Specify one of the following:"
303 : "\n\t- a single number"
304 : "\n\t- lists for blocks and/or block groups as follows:"
305 : "\n\t\t- [r, l_max] for spherical shell blocks and groups"
306 : "\n\t\t- [r, z] for filled cylinder blocks and groups containing them, "
307 : "\n\t\t where r must be > 2"
308 : "\n\t\t- [r, theta, z] for hollow cylinder blocks, where theta must be "
309 : "\n\t\t odd\n\n"
310 : "While the most verbose, the best choice for a production run is "
311 : "likely to specify a list for each cylindrical block instead of each "
312 : "cylindrical block group. If you set a whole group (e.g. InnerA) using "
313 : "[r, z], the interfaces between its filled cylinders "
314 : "(e.g. EAFilledCylinder) and its hollow cylinders (e.g. EACylinder) "
315 : "may not have similar resolution on each side unless r and z are "
316 : "close. This is because at these interfaces, the z direction in filled "
317 : "cylinders lines up with the radial direction in hollow cylinders. To "
318 : "get the resolution on either side of these interfaces to match well, "
319 : "you either want to set a cylindrical block group with r and z close "
320 : "in value or set the individual cylindrical blocks for more freedom. "
321 : "Also note that any h refinement in these cylindrical blocks will also "
322 : "be in similarly different directions at the interface, which affects "
323 : "this picture of trying to match the p refinement at the interface of "
324 : "hollow and filled cylinders."};
325 : };
326 :
327 0 : struct BoundaryConditions {
328 0 : static constexpr Options::String help = "The boundary conditions to apply.";
329 : };
330 : template <typename BoundaryConditionsBase>
331 0 : struct InnerBoundaryCondition {
332 0 : static std::string name() { return "InnerBoundary"; }
333 0 : static constexpr Options::String help =
334 : "Options for the inner boundary conditions.";
335 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
336 0 : using group = BoundaryConditions;
337 : };
338 :
339 : template <typename BoundaryConditionsBase>
340 0 : struct OuterBoundaryCondition {
341 0 : static std::string name() { return "OuterBoundary"; }
342 0 : static constexpr Options::String help =
343 : "Options for the outer boundary conditions.";
344 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
345 0 : using group = BoundaryConditions;
346 : };
347 :
348 0 : struct TimeDependentMaps {
349 0 : using type = Options::Auto<bco::TimeDependentMapOptions<true>,
350 : Options::AutoLabel::None>;
351 0 : static constexpr Options::String help =
352 : bco::TimeDependentMapOptions<true>::help;
353 : };
354 :
355 : template <typename Metavariables>
356 0 : using options = tmpl::append<
357 : tmpl::list<CenterA, CenterB, RadiusA, RadiusB, InitialRefinement,
358 : InitialGridPoints, OuterSphere, InnerSphereA, InnerSphereB,
359 : TimeDependentMaps>,
360 : tmpl::conditional_t<
361 : domain::BoundaryConditions::has_boundary_conditions_base_v<
362 : typename Metavariables::system>,
363 : tmpl::list<
364 : InnerBoundaryCondition<
365 : domain::BoundaryConditions::get_boundary_conditions_base<
366 : typename Metavariables::system>>,
367 : OuterBoundaryCondition<
368 : domain::BoundaryConditions::get_boundary_conditions_base<
369 : typename Metavariables::system>>>,
370 : tmpl::list<>>>;
371 :
372 0 : static constexpr Options::String help{
373 : "The CylindricalBinaryCompactObject domain is a general domain for "
374 : "two compact objects. The user must provide the (grid-frame) "
375 : "centers and radii of the spherical inner edge of the grid surrounding "
376 : "each of the two compact objects A and B."};
377 :
378 0 : CylindricalBinaryCompactObject(
379 : std::array<double, 3> center_A, std::array<double, 3> center_B,
380 : double radius_A, double radius_B,
381 : const typename InitialRefinement::type& initial_refinement,
382 : const typename InitialGridPoints::type& initial_grid_points,
383 : typename OuterSphere::type outer_sphere,
384 : std::optional<InnerSphere> inner_sphere_A = std::nullopt,
385 : std::optional<InnerSphere> inner_sphere_B = std::nullopt,
386 : std::optional<bco::TimeDependentMapOptions<true>> time_dependent_options =
387 : std::nullopt,
388 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
389 : inner_boundary_condition = nullptr,
390 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
391 : outer_boundary_condition = nullptr,
392 : const Options::Context& context = {});
393 :
394 0 : CylindricalBinaryCompactObject() = default;
395 0 : CylindricalBinaryCompactObject(const CylindricalBinaryCompactObject&) =
396 : delete;
397 0 : CylindricalBinaryCompactObject(CylindricalBinaryCompactObject&&) = default;
398 0 : CylindricalBinaryCompactObject& operator=(
399 : const CylindricalBinaryCompactObject&) = delete;
400 0 : CylindricalBinaryCompactObject& operator=(CylindricalBinaryCompactObject&&) =
401 : default;
402 0 : ~CylindricalBinaryCompactObject() override = default;
403 :
404 0 : Domain<3> create_domain() const override;
405 :
406 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
407 1 : grid_anchors() const override {
408 : return grid_anchors_;
409 : }
410 :
411 : std::vector<DirectionMap<
412 : 3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
413 1 : external_boundary_conditions() const override;
414 :
415 1 : std::vector<std::array<size_t, 3>> initial_extents() const override;
416 :
417 1 : std::vector<std::array<size_t, 3>> initial_refinement_levels() const override;
418 :
419 1 : auto functions_of_time(const std::unordered_map<std::string, double>&
420 : initial_expiration_times = {}) const
421 : -> std::unordered_map<
422 : std::string,
423 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>> override;
424 :
425 1 : std::vector<std::string> block_names() const override { return block_names_; }
426 :
427 : std::unordered_map<std::string, std::unordered_set<std::string>>
428 1 : block_groups() const override {
429 : return block_groups_;
430 : }
431 :
432 : private:
433 : // Note that center_A_ and center_B_ are rotated with respect to the
434 : // input centers (which are in the grid frame), so that we can
435 : // construct the map in a frame where the centers are offset in the
436 : // z direction. At the end, there will be another rotation back to
437 : // the grid frame (where the centers are offset in the x direction).
438 0 : std::array<double, 3> center_A_{};
439 0 : std::array<double, 3> center_B_{};
440 0 : double radius_A_{};
441 0 : double radius_B_{};
442 0 : double outer_radius_A_{};
443 0 : double outer_radius_B_{};
444 0 : typename std::vector<std::array<size_t, 3>> initial_refinement_{};
445 0 : typename std::vector<std::array<size_t, 3>> initial_grid_points_{};
446 0 : typename OuterSphere::type outer_sphere_{};
447 0 : std::optional<InnerSphere> inner_sphere_A_{};
448 0 : std::optional<InnerSphere> inner_sphere_B_{};
449 : // cut_spheres_offset_factor_ is eta in Eq. (A.9) of
450 : // https://arxiv.org/abs/1206.3015. cut_spheres_offset_factor_
451 : // could be set to unity to simplify the equations. Here we fix it
452 : // to the value 0.99 used in SpEC, so that we reproduce SpEC's
453 : // domain decomposition.
454 0 : double cut_spheres_offset_factor_{0.99};
455 : // z_cutting_plane_ is x_C in Eq. (A.9) of
456 : // https://arxiv.org/abs/1206.3015 (but rotated to the z-axis).
457 0 : double z_cutting_plane_{};
458 0 : size_t number_of_blocks_{};
459 0 : std::unordered_map<std::string, size_t> block_positions_{};
460 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
461 0 : inner_boundary_condition_;
462 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
463 0 : outer_boundary_condition_;
464 0 : std::vector<std::string> block_names_{};
465 : std::unordered_map<std::string, std::unordered_set<std::string>>
466 0 : block_groups_{};
467 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
468 0 : grid_anchors_{};
469 : // FunctionsOfTime options
470 0 : std::optional<bco::TimeDependentMapOptions<true>> time_dependent_options_{};
471 : };
472 : } // namespace domain::creators
|
