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 : class Interval;
36 : template <typename Map1, typename Map2>
37 : class ProductOf2Maps;
38 : template <typename Map1, typename Map2, typename Map3>
39 : class ProductOf3Maps;
40 : class SphericalToCartesianPfaffian;
41 : template <size_t VolumeDim>
42 : class Wedge;
43 : template <size_t VolumeDim>
44 : class DiscreteRotation;
45 : class UniformCylindricalEndcap;
46 : class UniformCylindricalFlatEndcap;
47 : class UniformCylindricalSide;
48 : } // namespace CoordinateMaps
49 :
50 : template <typename SourceFrame, typename TargetFrame, typename... Maps>
51 : class CoordinateMap;
52 :
53 : template <typename T>
54 : struct ExpandOverBlocks;
55 :
56 : namespace FunctionsOfTime {
57 : class FunctionOfTime;
58 : } // namespace FunctionsOfTime
59 : } // namespace domain
60 :
61 : namespace Frame {
62 : struct Grid;
63 : struct Distorted;
64 : struct Inertial;
65 : struct BlockLogical;
66 : } // namespace Frame
67 : /// \endcond
68 :
69 : namespace domain::creators {
70 :
71 : /*!
72 : * \ingroup ComputationalDomainGroup
73 : *
74 : * \brief A general domain for two compact objects based on cylinders.
75 : *
76 : * Creates a 3D Domain that represents a binary compact object
77 : * solution. This domain is described briefly in the Appendix of
78 : * \cite Buchman:2012dw, and is illustrated in Figure 20 of that
79 : * paper.
80 : *
81 : * In the code and options below, `ObjectA` and `ObjectB` refer to the
82 : * two compact objects. In the grid frame, `ObjectA` is located to the
83 : * right of (i.e. a more positive value of the x-coordinate than)
84 : * `ObjectB`. The inner edge of the Blocks surrounding each of
85 : * `ObjectA` and `ObjectB` is spherical in grid coordinates; the
86 : * user must specify the center and radius of this surface for both
87 : * `ObjectA` and `ObjectB`, and the user must specify the outer boundary
88 : * radius. The outer boundary is a sphere centered at the origin.
89 : *
90 : * This domain offers some grid anchors. See
91 : * `domain::creators::bco::create_grid_anchors` for which ones are offered.
92 : *
93 : * Note that Figure 20 of \cite Buchman:2012dw illustrates additional
94 : * spherical shells inside the "EA" and "EB" blocks, and the caption
95 : * of Figure 20 indicates that there are additional spherical shells
96 : * outside the "CA" and "CB" blocks; `CylindricalBinaryCompactObject`
97 : * has these extra shells inside "EA" only if the option `IncludeInnerSphereA`
98 : * is true, and it has the extra shells inside "EB" only if the option
99 : * `IncludeInnerSphereB` is true. If the shells are absent, then the "EA" and
100 : * "EB" blocks extend to the excision boundaries.
101 : *
102 : * The Blocks are named as follows:
103 : * - Each of CAFilledCylinder, EAFilledCylinder, EBFilledCylinder,
104 : * MAFilledCylinder, MBFilledCylinder, and CBFilledCylinder consists
105 : * of 5 blocks, named 'Center', 'East', 'North', 'West', and
106 : * 'South', so an example of a valid block name is
107 : * 'CAFilledCylinderCenter'.
108 : * - Each of CACylinder, EACylinder, EBCylinder, and CBCylinder
109 : * consists of 4 blocks, named 'East', 'North', 'West', and 'South',
110 : * so an example of a valid block name is 'CACylinderEast'.
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 maps_list = tmpl::flatten<
153 : tmpl::list<domain::CoordinateMap<
154 : Frame::BlockLogical, Frame::Inertial,
155 : CoordinateMaps::ProductOf3Maps<CoordinateMaps::Interval,
156 : CoordinateMaps::Interval,
157 : CoordinateMaps::Interval>,
158 : CoordinateMaps::UniformCylindricalEndcap,
159 : CoordinateMaps::DiscreteRotation<3>>,
160 : domain::CoordinateMap<
161 : Frame::BlockLogical, Frame::Inertial,
162 : CoordinateMaps::ProductOf2Maps<CoordinateMaps::Wedge<2>,
163 : CoordinateMaps::Interval>,
164 : CoordinateMaps::UniformCylindricalEndcap,
165 : CoordinateMaps::DiscreteRotation<3>>,
166 : domain::CoordinateMap<
167 : Frame::BlockLogical, Frame::Inertial,
168 : CoordinateMaps::ProductOf3Maps<CoordinateMaps::Interval,
169 : CoordinateMaps::Interval,
170 : CoordinateMaps::Interval>,
171 : CoordinateMaps::UniformCylindricalFlatEndcap,
172 : CoordinateMaps::DiscreteRotation<3>>,
173 : domain::CoordinateMap<
174 : Frame::BlockLogical, Frame::Inertial,
175 : CoordinateMaps::ProductOf2Maps<CoordinateMaps::Wedge<2>,
176 : CoordinateMaps::Interval>,
177 : CoordinateMaps::UniformCylindricalFlatEndcap,
178 : CoordinateMaps::DiscreteRotation<3>>,
179 : domain::CoordinateMap<
180 : Frame::BlockLogical, Frame::Inertial,
181 : CoordinateMaps::ProductOf2Maps<CoordinateMaps::Wedge<2>,
182 : CoordinateMaps::Interval>,
183 : CoordinateMaps::UniformCylindricalSide,
184 : CoordinateMaps::DiscreteRotation<3>>,
185 : domain::CoordinateMap<
186 : Frame::BlockLogical, Frame::Inertial,
187 : domain::CoordinateMaps::ProductOf2Maps<
188 : CoordinateMaps::Interval, CoordinateMaps::Identity<2>>,
189 : domain::CoordinateMaps::SphericalToCartesianPfaffian,
190 : CoordinateMaps::ProductOf3Maps<CoordinateMaps::Affine,
191 : CoordinateMaps::Affine,
192 : CoordinateMaps::Affine>>,
193 : bco::TimeDependentMapOptions<true>::maps_list>>;
194 :
195 0 : struct CenterA {
196 0 : using type = std::array<double, 3>;
197 0 : static constexpr Options::String help = {
198 : "Grid coordinates of center for Object A, which is at x>0."};
199 : };
200 0 : struct CenterB {
201 0 : using type = std::array<double, 3>;
202 0 : static constexpr Options::String help = {
203 : "Grid coordinates of center for Object B, which is at x<0."};
204 : };
205 0 : struct RadiusA {
206 0 : using type = double;
207 0 : static constexpr Options::String help = {
208 : "Grid-coordinate radius of grid boundary around Object A."};
209 : };
210 0 : struct RadiusB {
211 0 : using type = double;
212 0 : static constexpr Options::String help = {
213 : "Grid-coordinate radius of grid boundary around Object B."};
214 : };
215 0 : struct IncludeInnerSphereA {
216 0 : using type = bool;
217 0 : static constexpr Options::String help = {
218 : "Add an extra spherical layer of Blocks around Object A."};
219 : };
220 0 : struct IncludeInnerSphereB {
221 0 : using type = bool;
222 0 : static constexpr Options::String help = {
223 : "Add an extra spherical layer of Blocks around Object B."};
224 : };
225 0 : struct OuterRadius {
226 0 : using type = double;
227 0 : static constexpr Options::String help = {
228 : "Grid-coordinate radius of outer boundary."};
229 : };
230 0 : struct UseEquiangularMap {
231 0 : using type = bool;
232 0 : static constexpr Options::String help = {
233 : "Distribute grid points equiangularly in 2d wedges."};
234 0 : static bool suggested_value() { return false; }
235 : };
236 :
237 0 : struct InitialRefinement {
238 0 : using type = std::variant<
239 : size_t, std::array<size_t, 3>, std::vector<std::array<size_t, 3>>,
240 : std::unordered_map<std::string, std::array<size_t, 3>>,
241 : std::unordered_map<std::string,
242 : std::variant<std::array<size_t, 3>, size_t>>>;
243 0 : static constexpr Options::String help = {
244 : "Initial refinement level. Specify one of: a single number, a list "
245 : "representing [r, theta, perp], or such a list for every block in the "
246 : "domain. Here 'r' is the radial direction normal to the inner and "
247 : "outer boundaries, 'theta' is the periodic direction, and 'perp' is "
248 : "the third direction. Note that for spherical shell block groups "
249 : "('InnerSphereA', 'InnerSphereB', and 'OuterSphere'), you must instead "
250 : "specify refinement as a single value representing radial refinement."};
251 : };
252 0 : struct InitialGridPoints {
253 0 : using type = std::variant<
254 : size_t, std::array<size_t, 3>, std::vector<std::array<size_t, 3>>,
255 : std::unordered_map<std::string, std::array<size_t, 3>>,
256 : std::unordered_map<std::string, std::variant<std::array<size_t, 3>,
257 : std::array<size_t, 2>>>>;
258 0 : static constexpr Options::String help = {
259 : "Initial number of grid points. Specify one of: a single number, a "
260 : "list representing [r, theta, perp], or such a list for every block in "
261 : "the domain. Here 'r' is the radial direction normal to the inner and "
262 : "outer boundaries, 'theta' is the periodic direction, and 'perp' is "
263 : "the third direction. The exception to this is that for spherical "
264 : "shell blocks groups ('InnerSphereA', 'InnerSphereB', 'OuterSphere'),"
265 : "you must instead specify grid points as [r, L_max]. The exception to "
266 : "this is if a single number is specified for global initial grid "
267 : "points."};
268 : };
269 :
270 0 : struct BoundaryConditions {
271 0 : static constexpr Options::String help = "The boundary conditions to apply.";
272 : };
273 : template <typename BoundaryConditionsBase>
274 0 : struct InnerBoundaryCondition {
275 0 : static std::string name() { return "InnerBoundary"; }
276 0 : static constexpr Options::String help =
277 : "Options for the inner boundary conditions.";
278 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
279 0 : using group = BoundaryConditions;
280 : };
281 :
282 : template <typename BoundaryConditionsBase>
283 0 : struct OuterBoundaryCondition {
284 0 : static std::string name() { return "OuterBoundary"; }
285 0 : static constexpr Options::String help =
286 : "Options for the outer boundary conditions.";
287 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
288 0 : using group = BoundaryConditions;
289 : };
290 :
291 0 : struct TimeDependentMaps {
292 0 : using type = Options::Auto<bco::TimeDependentMapOptions<true>,
293 : Options::AutoLabel::None>;
294 0 : static constexpr Options::String help =
295 : bco::TimeDependentMapOptions<true>::help;
296 : };
297 :
298 : template <typename Metavariables>
299 0 : using options = tmpl::append<
300 : tmpl::list<CenterA, CenterB, RadiusA, RadiusB, IncludeInnerSphereA,
301 : IncludeInnerSphereB, OuterRadius, UseEquiangularMap,
302 : InitialRefinement, InitialGridPoints, TimeDependentMaps>,
303 : tmpl::conditional_t<
304 : domain::BoundaryConditions::has_boundary_conditions_base_v<
305 : typename Metavariables::system>,
306 : tmpl::list<
307 : InnerBoundaryCondition<
308 : domain::BoundaryConditions::get_boundary_conditions_base<
309 : typename Metavariables::system>>,
310 : OuterBoundaryCondition<
311 : domain::BoundaryConditions::get_boundary_conditions_base<
312 : typename Metavariables::system>>>,
313 : tmpl::list<>>>;
314 :
315 0 : static constexpr Options::String help{
316 : "The CylindricalBinaryCompactObject domain is a general domain for "
317 : "two compact objects. The user must provide the (grid-frame) "
318 : "centers and radii of the spherical inner edge of the grid surrounding "
319 : "each of the two compact objects A and B."};
320 :
321 0 : CylindricalBinaryCompactObject(
322 : std::array<double, 3> center_A, std::array<double, 3> center_B,
323 : double radius_A, double radius_B, bool include_inner_sphere_A,
324 : bool include_inner_sphere_B, double outer_radius,
325 : bool use_equiangular_map,
326 : const typename InitialRefinement::type& initial_refinement,
327 : const typename InitialGridPoints::type& initial_grid_points,
328 : std::optional<bco::TimeDependentMapOptions<true>> time_dependent_options =
329 : std::nullopt,
330 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
331 : inner_boundary_condition = nullptr,
332 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
333 : outer_boundary_condition = nullptr,
334 : const Options::Context& context = {});
335 :
336 0 : CylindricalBinaryCompactObject() = default;
337 0 : CylindricalBinaryCompactObject(const CylindricalBinaryCompactObject&) =
338 : delete;
339 0 : CylindricalBinaryCompactObject(CylindricalBinaryCompactObject&&) = default;
340 0 : CylindricalBinaryCompactObject& operator=(
341 : const CylindricalBinaryCompactObject&) = delete;
342 0 : CylindricalBinaryCompactObject& operator=(CylindricalBinaryCompactObject&&) =
343 : default;
344 0 : ~CylindricalBinaryCompactObject() override = default;
345 :
346 0 : Domain<3> create_domain() const override;
347 :
348 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
349 1 : grid_anchors() const override {
350 : return grid_anchors_;
351 : }
352 :
353 : std::vector<DirectionMap<
354 : 3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
355 1 : external_boundary_conditions() const override;
356 :
357 1 : std::vector<std::array<size_t, 3>> initial_extents() const override;
358 :
359 1 : std::vector<std::array<size_t, 3>> initial_refinement_levels() const override;
360 :
361 1 : auto functions_of_time(const std::unordered_map<std::string, double>&
362 : initial_expiration_times = {}) const
363 : -> std::unordered_map<
364 : std::string,
365 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>> override;
366 :
367 1 : std::vector<std::string> block_names() const override { return block_names_; }
368 :
369 : std::unordered_map<std::string, std::unordered_set<std::string>>
370 1 : block_groups() const override {
371 : return block_groups_;
372 : }
373 :
374 : private:
375 : // Note that center_A_ and center_B_ are rotated with respect to the
376 : // input centers (which are in the grid frame), so that we can
377 : // construct the map in a frame where the centers are offset in the
378 : // z direction. At the end, there will be another rotation back to
379 : // the grid frame (where the centers are offset in the x direction).
380 0 : std::array<double, 3> center_A_{};
381 0 : std::array<double, 3> center_B_{};
382 0 : double radius_A_{};
383 0 : double radius_B_{};
384 0 : double outer_radius_A_{};
385 0 : double outer_radius_B_{};
386 0 : bool include_inner_sphere_A_{};
387 0 : bool include_inner_sphere_B_{};
388 0 : double outer_radius_{};
389 0 : bool use_equiangular_map_{false};
390 0 : typename std::vector<std::array<size_t, 3>> initial_refinement_{};
391 0 : typename std::vector<std::array<size_t, 3>> initial_grid_points_{};
392 : // cut_spheres_offset_factor_ is eta in Eq. (A.9) of
393 : // https://arxiv.org/abs/1206.3015. cut_spheres_offset_factor_
394 : // could be set to unity to simplify the equations. Here we fix it
395 : // to the value 0.99 used in SpEC, so that we reproduce SpEC's
396 : // domain decomposition.
397 0 : double cut_spheres_offset_factor_{0.99};
398 : // z_cutting_plane_ is x_C in Eq. (A.9) of
399 : // https://arxiv.org/abs/1206.3015 (but rotated to the z-axis).
400 0 : double z_cutting_plane_{};
401 0 : size_t number_of_blocks_{};
402 0 : size_t first_outer_shell_block{};
403 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
404 0 : inner_boundary_condition_;
405 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
406 0 : outer_boundary_condition_;
407 0 : std::vector<std::string> block_names_{};
408 : std::unordered_map<std::string, std::unordered_set<std::string>>
409 0 : block_groups_{};
410 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
411 0 : grid_anchors_{};
412 : // FunctionsOfTime options
413 0 : std::optional<bco::TimeDependentMapOptions<true>> time_dependent_options_{};
414 : };
415 : } // namespace domain::creators
|
