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