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