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 <limits>
9 : #include <memory>
10 : #include <optional>
11 : #include <string>
12 : #include <type_traits>
13 : #include <unordered_map>
14 : #include <unordered_set>
15 : #include <variant>
16 : #include <vector>
17 :
18 : #include "DataStructures/Tensor/Tensor.hpp"
19 : #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
20 : #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp"
21 : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
22 : #include "Domain/CoordinateMaps/Distribution.hpp"
23 : #include "Domain/Creators/DomainCreator.hpp"
24 : #include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp"
25 : #include "Domain/Domain.hpp"
26 : #include "Domain/Structure/DirectionMap.hpp"
27 : #include "Options/Auto.hpp"
28 : #include "Options/Context.hpp"
29 : #include "Options/String.hpp"
30 : #include "Utilities/ErrorHandling/Assert.hpp"
31 : #include "Utilities/TMPL.hpp"
32 :
33 : /// \cond
34 : namespace domain {
35 : namespace CoordinateMaps {
36 : class Affine;
37 : class Equiangular;
38 : template <size_t VolumeDim>
39 : class Identity;
40 : class Interval;
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 Dim>
47 : class Wedge;
48 : template <size_t VolumeDim>
49 : class DiscreteRotation;
50 : class Frustum;
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 : * \ingroup ComputationalDomainGroup
75 : *
76 : * \brief A general domain for two compact objects.
77 : *
78 : * \image html binary_compact_object_domain.png "A BHNS domain."
79 : *
80 : * Creates a 3D Domain that represents a binary compact object solution. The
81 : * Domain consists of 4 or 5 nested layers of blocks; these layers are, working
82 : * from the interior toward the exterior:
83 : *
84 : * - **Object A/B interior**: (optional) The block at the center of each
85 : * compact object, if not excised. If present, this block is a cube. If
86 : * excised, the hole left by its absence is spherical.
87 : * - **Object A/B shell**: The 6 blocks that resolve each individual compact
88 : * object. This layer has a spherical outer boundary - if the corresponding
89 : * interior block exists, then the layer is a cube-to-sphere transition; if
90 : * the interior block is excised, then the layer is a spherical shell.
91 : * - **Object A/B cube**: The 6 blocks that surround each object with a cube.
92 : * Around each compact object, this layer transitions from a sphere to a cube.
93 : * - **Envelope**: The 10 blocks that transition from the two inner cubes to a
94 : * sphere centered at the origin.
95 : * - **Outer shell**: The 10 blocks that form an outer shell centered at the
96 : * origin, consisting of 2 endcap Wedges on the +x and -x axes, and 8 half
97 : * Wedges along the yz plane. This layer is spherical, so a logarithmic map
98 : * can optionally be used in this layer. This allows the domain to extend to
99 : * large radial distances from the compact objects. This layer can be
100 : * h-refined radially, creating a layer of multiple concentric spherical
101 : * shells.
102 : *
103 : * \par Notes:
104 : * - Object A is located to the right of the origin (along the positive x-axis)
105 : * and Object B is located to the left of the origin in the Grid frame.
106 : * - This domain offers some grid anchors. See
107 : * `domain::creators::bco::create_grid_anchors` for which ones are offered.
108 : * - "Cutting plane" refers to the plane along which the domain divides into two
109 : * hemispheres. The cutting plane always intersects the x-axis at the origin.
110 : * - The x-coordinate locations of the two objects should be chosen such that
111 : * the center of mass is located at x=0 at the initial time (typically t=0).
112 : * - The cubes are first constructed at the origin. Then, they are translated
113 : * left/right by their Object's x-coordinate and offset depending on the cube
114 : * length.
115 : * - The CubeScale option describes how to scale the length of the cube
116 : * surrounding object A/B. It must be greater than or equal to 1.0 with 1.0
117 : * meaning the side length of the cube is the initial physical separation
118 : * between the two objects. If CubeScale is greater than 1.0, the centers of
119 : * the two objects will be offset relative to the centers of the cubes.
120 : * - Alternatively, one can replace the inner shell and cube blocks of each
121 : * object with a single cartesian cube. This is less efficient, but allows
122 : * testing of methods only coded on cartesian grids.
123 : *
124 : * \par Time dependence:
125 : * The following time-dependent maps are applied:
126 : *
127 : * - A piecewise `Expansion`, a `Rotation` and a piecewise `Translation` is
128 : * applied to all blocks from the Grid to the Inertial frame. However, if there
129 : * is a shape map in the block (defined below), then the expansion, rotation,
130 : * and translation maps go from the Distorted to the Inertial frame.
131 : * - If an object is excised, then the corresponding shell has a
132 : * `Shape` map. The shape map goes from the Grid to the Distorted frame.
133 : *
134 : * All time dependent maps are optional to specify. To include a map, specify
135 : * its options. Otherwise specify `None` for that map. You can also turn off
136 : * time dependent maps all together by specifying `None` for the
137 : * `TimeDependentMaps` option. See
138 : * `domain::creators::bco::TimeDependentMapOptions`. This class must pass a
139 : * template parameter of `false` to
140 : * `domain::creators::bco::TimeDependentMapOptions`.
141 : *
142 : * The wavezone can support using spherical harmonic shells by setting
143 : * `SphericalHarmonicsInWavezone` to `true`. In this case you cannot set an
144 : * angular h-refinement level in the angular directions, and you can only choose
145 : * the \f$\ell\f$ of the shells since having an \f$m_{\max}\f$ below the maximum
146 : * allowed from the \f$\ell\f$ can drive simulations unstable.
147 : *
148 : * The `UseWorldtube` option defaults to `false`. When set to `true`, some of
149 : * the functions of time will be `IntegratedFunctionOfTime` used to control the
150 : * orbit of the worldtube.
151 : */
152 1 : class BinaryCompactObject : public DomainCreator<3> {
153 : private:
154 : // Time-independent maps
155 0 : using Affine = CoordinateMaps::Affine;
156 0 : using Affine3D = CoordinateMaps::ProductOf3Maps<Affine, Affine, Affine>;
157 0 : using Identity2D = CoordinateMaps::Identity<2>;
158 : // The Translation type is no longer needed, but it is kept here for backwards
159 : // compatibility with old domains.
160 0 : using Translation = CoordinateMaps::ProductOf2Maps<Affine, Identity2D>;
161 0 : using Equiangular = CoordinateMaps::Equiangular;
162 0 : using Equiangular3D =
163 : CoordinateMaps::ProductOf3Maps<Equiangular, Equiangular, Equiangular>;
164 :
165 : public:
166 0 : using maps_list = tmpl::flatten<tmpl::list<
167 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Affine3D>,
168 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
169 : Equiangular3D>,
170 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Affine3D,
171 : Translation>,
172 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
173 : CoordinateMaps::DiscreteRotation<3>, Affine3D>,
174 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
175 : Equiangular3D>,
176 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Equiangular3D,
177 : Translation>,
178 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
179 : CoordinateMaps::Frustum>,
180 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
181 : CoordinateMaps::Wedge<3>>,
182 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
183 : CoordinateMaps::Wedge<3>, Translation>,
184 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Affine3D,
185 : Affine3D>,
186 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Equiangular3D,
187 : Affine3D>,
188 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
189 : CoordinateMaps::Wedge<3>, Affine3D>,
190 : domain::CoordinateMap<
191 : Frame::BlockLogical, Frame::Inertial,
192 : domain::CoordinateMaps::ProductOf2Maps<
193 : domain::CoordinateMaps::Interval,
194 : domain::CoordinateMaps::Identity<2>>,
195 : domain::CoordinateMaps::SphericalToCartesianPfaffian>,
196 : bco::TimeDependentMapOptions<false>::maps_list>>;
197 :
198 : /// Options for an excision region in the domain
199 1 : struct Excision {
200 0 : static constexpr Options::String help = {
201 : "Excise the interior of the object, leaving a spherical hole in its "
202 : "absence."};
203 : template <typename BoundaryConditionsBase>
204 0 : struct BoundaryCondition {
205 0 : static std::string name() { return "ExciseWithBoundaryCondition"; }
206 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
207 0 : static constexpr Options::String help = {
208 : "The boundary condition to impose on the excision surface."};
209 : };
210 : template <typename Metavariables>
211 0 : using options = tmpl::list<BoundaryCondition<
212 : domain::BoundaryConditions::get_boundary_conditions_base<
213 : typename Metavariables::system>>>;
214 0 : Excision() = default;
215 : // NOLINTNEXTLINE(google-explicit-constructor)
216 0 : Excision(std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
217 : boundary_condition_in)
218 : : boundary_condition(std::move(boundary_condition_in)) {}
219 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
220 0 : boundary_condition;
221 : };
222 :
223 : /// Options for one of the two objects in the binary domain
224 1 : struct Object {
225 0 : static constexpr Options::String help = {
226 : "Options for an object in a binary domain."};
227 0 : struct InnerRadius {
228 0 : using type = double;
229 0 : static constexpr Options::String help = {
230 : "Inner coordinate radius of Layer 1."};
231 0 : static double lower_bound() { return 0.; }
232 : };
233 0 : struct OuterRadius {
234 0 : using type = double;
235 0 : static constexpr Options::String help = {
236 : "Outer coordinate radius of Layer 1"};
237 0 : static double lower_bound() { return 0.; }
238 : };
239 0 : struct XCoord {
240 0 : using type = double;
241 0 : static constexpr Options::String help = {"x-coordinate of center."};
242 : };
243 0 : struct Interior {
244 0 : using type = Options::Auto<Excision>;
245 0 : static constexpr Options::String help = {
246 : "Specify 'ExciseWithBoundaryCondition' and a boundary condition to "
247 : "excise Layer 0, leaving a spherical hole in its absence, or set to "
248 : "'Auto' to fill the interior."};
249 : };
250 0 : struct ExciseInterior {
251 0 : using type = bool;
252 0 : static constexpr Options::String help = {
253 : "Excise Layer 0, leaving a spherical hole in its absence."};
254 : };
255 0 : struct UseLogarithmicMap {
256 0 : using type = bool;
257 0 : static constexpr Options::String help = {
258 : "Use a logarithmically spaced radial grid in the part of Layer 1 "
259 : "enveloping the object (requires the interior is excised)"};
260 : };
261 : template <typename Metavariables>
262 0 : using options = tmpl::list<
263 : InnerRadius, OuterRadius, XCoord,
264 : tmpl::conditional_t<
265 : domain::BoundaryConditions::has_boundary_conditions_base_v<
266 : typename Metavariables::system>,
267 : Interior, ExciseInterior>,
268 : UseLogarithmicMap>;
269 0 : Object() {} // NOLINT(modernize-use-equals-default)
270 0 : Object(double local_inner_radius, double local_outer_radius,
271 : double local_x_coord, std::optional<Excision> interior,
272 : bool local_use_logarithmic_map)
273 : : inner_radius(local_inner_radius),
274 : outer_radius(local_outer_radius),
275 : x_coord(local_x_coord),
276 : inner_boundary_condition(
277 : interior.has_value()
278 : ? std::make_optional(std::move(interior->boundary_condition))
279 : : std::nullopt),
280 : use_logarithmic_map(local_use_logarithmic_map) {}
281 0 : Object(double local_inner_radius, double local_outer_radius,
282 : double local_x_coord, bool local_excise_interior,
283 : bool local_use_logarithmic_map)
284 : : inner_radius(local_inner_radius),
285 : outer_radius(local_outer_radius),
286 : x_coord(local_x_coord),
287 : inner_boundary_condition(
288 : local_excise_interior
289 : ? std::optional<std::unique_ptr<
290 : domain::BoundaryConditions::BoundaryCondition>>{nullptr}
291 : : std::nullopt),
292 : use_logarithmic_map(local_use_logarithmic_map) {}
293 :
294 : /// Whether or not the object should be excised from the domain, leaving a
295 : /// spherical hole. When this is true, `inner_boundary_condition` is
296 : /// guaranteed to hold a value (though it might be a `nullptr` if we are not
297 : /// working with boundary conditions).
298 1 : bool is_excised() const;
299 :
300 0 : double inner_radius{};
301 0 : double outer_radius{};
302 0 : double x_coord{};
303 : std::optional<
304 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>
305 0 : inner_boundary_condition;
306 0 : bool use_logarithmic_map{};
307 : };
308 :
309 : // Simpler version of an object: a single cube centered on (xCoord,0,0)
310 0 : struct CartesianCubeAtXCoord {
311 0 : static constexpr Options::String help = {
312 : "Options to set a single cube at a location on the x-axis"};
313 0 : struct XCoord {
314 0 : static std::string name() { return "CartesianCubeAtXCoord"; }
315 0 : using type = double;
316 0 : static constexpr Options::String help = {"x-coordinate of center."};
317 : };
318 0 : using options = tmpl::list<XCoord>;
319 0 : CartesianCubeAtXCoord() = default;
320 : // NOLINTNEXTLINE(google-explicit-constructor)
321 0 : CartesianCubeAtXCoord(const double x_coord_in) : x_coord(x_coord_in) {}
322 0 : bool is_excised() const { return false; }
323 0 : double x_coord;
324 : };
325 :
326 0 : struct ObjectA {
327 0 : using type = std::variant<Object, CartesianCubeAtXCoord>;
328 0 : static constexpr Options::String help = {
329 : "Options for the object to the right of the origin (along the positive "
330 : "x-axis)."};
331 : };
332 :
333 0 : struct ObjectB {
334 0 : using type = std::variant<Object, CartesianCubeAtXCoord>;
335 0 : static constexpr Options::String help = {
336 : "Options for the object to the left of the origin (along the negative "
337 : "x-axis)."};
338 : };
339 :
340 0 : struct CenterOfMassOffset {
341 0 : using type = std::array<double, 2>;
342 0 : static constexpr Options::String help = {
343 : "Offset in the y and z axes applied to both object A and B in order to "
344 : "control the center of mass. This moves the location of the two objects"
345 : " in the grid frame but keeps the Envelope and OuterShell centered on "
346 : "the origin in the grid frame."};
347 : };
348 :
349 0 : struct Envelope {
350 0 : static constexpr Options::String help = {
351 : "Options for the sphere enveloping the two objects."};
352 : };
353 :
354 0 : struct EnvelopeRadius {
355 0 : using group = Envelope;
356 0 : static std::string name() { return "Radius"; }
357 0 : using type = double;
358 0 : static constexpr Options::String help = {
359 : "Radius of the sphere enveloping the two objects."};
360 : };
361 :
362 0 : struct OuterShell {
363 0 : static constexpr Options::String help = {
364 : "Options for the outer spherical shell."};
365 : };
366 :
367 0 : struct OuterRadius {
368 0 : using group = OuterShell;
369 0 : static std::string name() { return "Radius"; }
370 0 : using type = double;
371 0 : static constexpr Options::String help = {"Radius of the entire domain."};
372 : };
373 :
374 0 : struct RadialPartitioningOuterShell {
375 0 : static std::string name() { return "RadialPartitioning"; }
376 0 : using group = OuterShell;
377 0 : using type = std::vector<double>;
378 0 : static constexpr Options::String help = {
379 : "Radial coordinates of the boundaries splitting the outer spherical "
380 : "shells between the envelope radius and OuterRadius. They must be "
381 : "given in ascending order. This should be used if boundaries need to "
382 : "be set at specific radii. If the number but not the specific "
383 : "locations of the boundaries are important, use InitialRefinement "
384 : "instead."};
385 : };
386 :
387 0 : struct RadialDistributionOuterShell {
388 0 : static std::string name() { return "RadialDistribution"; }
389 0 : using group = OuterShell;
390 0 : using type =
391 : std::variant<domain::CoordinateMaps::Distribution,
392 : std::vector<domain::CoordinateMaps::Distribution>>;
393 0 : static constexpr Options::String help = {
394 : "Select the radial distribution of grid points in each outer spherical "
395 : "shell. There must be N+1 radial distributions specified for N radial "
396 : "partitions. You can also specify just a single radial distribution "
397 : "(not in a vector) which will use the same distribution for all "
398 : "partitions."};
399 : };
400 :
401 0 : struct OpeningAngle {
402 0 : using group = OuterShell;
403 0 : static std::string name() { return "OpeningAngle"; }
404 0 : using type = double;
405 0 : static constexpr Options::String help = {
406 : "The combined opening angle of the two half wedges of the outer shell"
407 : " in degrees. A value of 120.0 partitions the x-y and x-z slices of the"
408 : " outer shell into six Blocks of equal angular size."};
409 : };
410 :
411 0 : struct SphericalHarmonicsInWavezone {
412 0 : using group = OuterShell;
413 0 : static std::string name() { return "UseSphericalHarmonics"; }
414 0 : using type = bool;
415 0 : static bool suggested_value() { return false; }
416 0 : static constexpr Options::String help = {
417 : "Use a spherical-harmonic basis for the outer wavezone shell(s) "
418 : "instead of the default 10-wedge Cartesian basis. When enabled, "
419 : "InitialGridPoints for wavezone blocks must be specified as "
420 : "array<size_t, 2> = {radial_points, L_max}, and InitialRefinement "
421 : "as a scalar or single-element value (angular refinement is fixed at "
422 : "0). "};
423 : };
424 :
425 0 : struct UseWorldtube {
426 0 : using type = bool;
427 0 : static constexpr Options::String help = {
428 : "Whether to set up functions of time appropriate for a worldtube run. "
429 : "When true, some functions of time will be IntegratedFunctionOfTime "
430 : "used to control the orbit of the worldtube. Used by the curved scalar "
431 : "wave worldtube executable."};
432 : };
433 :
434 0 : struct CubeScale {
435 0 : using type = double;
436 0 : static constexpr Options::String help = {
437 : "Specify the desired cube scale that must be greater than or equal to "
438 : "1.0. The initial separation is multiplied by this cube scale to "
439 : "produce larger cubes around each object which is desirable when "
440 : "closer to merger."};
441 0 : static double lower_bound() { return 1.0; }
442 : };
443 :
444 0 : struct InitialRefinement {
445 0 : using type = std::variant<
446 : size_t, std::array<size_t, 3>, std::vector<std::array<size_t, 3>>,
447 : std::unordered_map<std::string, std::array<size_t, 3>>,
448 : std::unordered_map<std::string,
449 : std::variant<std::array<size_t, 3>, size_t>>>;
450 0 : static constexpr Options::String help = {
451 : "Initial refinement level in each block of the domain. See main help "
452 : "text for details."};
453 : };
454 :
455 0 : struct InitialGridPoints {
456 0 : using type = std::variant<
457 : size_t, std::array<size_t, 3>, std::vector<std::array<size_t, 3>>,
458 : std::unordered_map<std::string, std::array<size_t, 3>>,
459 : std::unordered_map<std::string, std::variant<std::array<size_t, 3>,
460 : std::array<size_t, 2>>>>;
461 0 : static constexpr Options::String help = {
462 : "Initial number of grid points in the elements of each block of the "
463 : "domain. See main help text for details."};
464 : };
465 :
466 0 : struct UseEquiangularMap {
467 0 : using type = bool;
468 0 : static constexpr Options::String help = {
469 : "Distribute grid points equiangularly."};
470 0 : static bool suggested_value() { return true; }
471 : };
472 :
473 0 : struct RadialDistributionEnvelope {
474 0 : using group = Envelope;
475 0 : static std::string name() { return "RadialDistribution"; }
476 0 : using type = CoordinateMaps::Distribution;
477 0 : static constexpr Options::String help = {
478 : "The distribution of radial grid points in the envelope, the layer "
479 : "made of ten bulged Frustums."};
480 : };
481 :
482 : template <typename BoundaryConditionsBase>
483 0 : struct OuterBoundaryCondition {
484 0 : using group = OuterShell;
485 0 : static std::string name() { return "BoundaryCondition"; }
486 0 : static constexpr Options::String help =
487 : "Options for the outer boundary conditions.";
488 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
489 : };
490 :
491 : // This is for optional time dependent maps
492 0 : struct TimeDependentMaps {
493 0 : using type = Options::Auto<bco::TimeDependentMapOptions<false>,
494 : Options::AutoLabel::None>;
495 0 : static constexpr Options::String help =
496 : bco::TimeDependentMapOptions<false>::help;
497 : };
498 :
499 : template <typename Metavariables>
500 0 : using options = tmpl::append<
501 : tmpl::list<ObjectA, ObjectB, CenterOfMassOffset, EnvelopeRadius,
502 : OuterRadius, CubeScale, InitialRefinement, InitialGridPoints,
503 : UseEquiangularMap, RadialDistributionEnvelope,
504 : RadialPartitioningOuterShell, RadialDistributionOuterShell,
505 : OpeningAngle, SphericalHarmonicsInWavezone, UseWorldtube,
506 : TimeDependentMaps>,
507 : tmpl::conditional_t<
508 : domain::BoundaryConditions::has_boundary_conditions_base_v<
509 : typename Metavariables::system>,
510 : tmpl::list<OuterBoundaryCondition<
511 : domain::BoundaryConditions::get_boundary_conditions_base<
512 : typename Metavariables::system>>>,
513 : tmpl::list<>>>;
514 :
515 0 : static constexpr Options::String help{
516 : "A general domain for two compact objects. Each object is represented by "
517 : "a cube along the x-axis. Object A is located on the right and Object B "
518 : "is located on the left. Their locations should be chosen such that "
519 : "their center of mass is located at the origin."
520 : "The interior of each object can have a spherical excision to "
521 : "represent a black hole."
522 : "\n"
523 : "The two objects are enveloped by a sphere centered at the origin, "
524 : "and by an outer shell that can transition to large outer radii."
525 : "\n"
526 : "Both the InitialRefinement and the InitialGridPoints can be one of "
527 : "the following:\n"
528 : " - A single number: Uniform refinement in all blocks and "
529 : "dimensions\n"
530 : " - Three numbers: Refinement in [polar, azimuthal, radial] direction "
531 : "in all blocks\n"
532 : " - A map from block names or groups to three numbers: Per-block "
533 : "refinement in [polar, azimuthal, radial] direction\n"
534 : " - A list, with [polar, azimuthal, radial] refinement for each block\n"
535 : "\n"
536 : "The domain can rotate around the "
537 : "z-axis and expand/compress radially. The two objects can each have a "
538 : "spherical distortion (shape map)."};
539 :
540 0 : BinaryCompactObject(
541 : typename ObjectA::type object_A, typename ObjectB::type object_B,
542 : std::array<double, 2> center_of_mass_offset, double envelope_radius,
543 : double outer_radius, double cube_scale,
544 : const typename InitialRefinement::type& initial_refinement,
545 : const typename InitialGridPoints::type& initial_number_of_grid_points,
546 : bool use_equiangular_map = true,
547 : CoordinateMaps::Distribution radial_distribution_envelope =
548 : CoordinateMaps::Distribution::Projective,
549 : const std::vector<double>& radial_partitioning_outer_shell = {},
550 : const typename RadialDistributionOuterShell::type&
551 : radial_distribution_outer_shell =
552 : CoordinateMaps::Distribution::Linear,
553 : double opening_angle_in_degrees = 90.0,
554 : bool spherical_harmonics_in_wavezone = false, bool use_worldtube = false,
555 : std::optional<bco::TimeDependentMapOptions<false>>
556 : time_dependent_options = std::nullopt,
557 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
558 : outer_boundary_condition = nullptr,
559 : const Options::Context& context = {});
560 :
561 0 : BinaryCompactObject() = default;
562 0 : BinaryCompactObject(const BinaryCompactObject&) = delete;
563 0 : BinaryCompactObject(BinaryCompactObject&&) = default;
564 0 : BinaryCompactObject& operator=(const BinaryCompactObject&) = delete;
565 0 : BinaryCompactObject& operator=(BinaryCompactObject&&) = default;
566 0 : ~BinaryCompactObject() override = default;
567 :
568 0 : Domain<3> create_domain() const override;
569 :
570 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
571 1 : grid_anchors() const override {
572 : return grid_anchors_;
573 : }
574 :
575 : std::vector<DirectionMap<
576 : 3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
577 1 : external_boundary_conditions() const override;
578 :
579 1 : std::vector<std::array<size_t, 3>> initial_extents() const override;
580 :
581 1 : std::vector<std::array<size_t, 3>> initial_refinement_levels()
582 : const override {
583 : return initial_refinement_;
584 : }
585 :
586 1 : std::vector<std::string> block_names() const override { return block_names_; }
587 :
588 : std::unordered_map<std::string, std::unordered_set<std::string>>
589 1 : block_groups() const override {
590 : return block_groups_;
591 : }
592 :
593 1 : auto functions_of_time(const std::unordered_map<std::string, double>&
594 : initial_expiration_times = {}) const
595 : -> std::unordered_map<
596 : std::string,
597 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>> override;
598 :
599 : private:
600 0 : typename ObjectA::type object_A_{Object{}};
601 0 : typename ObjectB::type object_B_{Object{}};
602 0 : std::array<double, 2> center_of_mass_offset_{};
603 0 : double envelope_radius_ = std::numeric_limits<double>::signaling_NaN();
604 0 : double outer_radius_ = std::numeric_limits<double>::signaling_NaN();
605 0 : std::vector<std::array<size_t, 3>> initial_refinement_;
606 : // For most blocks this stores the number of grid points in each logical
607 : // direction. For spherical-harmonic shell blocksthe angular entries instead
608 : // store the spherical-harmonic degrees {n_radial, l_max, m_max}, because
609 : // ell is unambiguous whereas the implied number of collocation points
610 : // depends on the spectral implementation. The conversion of (l_max, m_max)
611 : // to the actual number of collocation points is applied in
612 : // `initial_extents()`.
613 0 : std::vector<std::array<size_t, 3>> initial_number_of_grid_points_;
614 0 : bool use_equiangular_map_ = true;
615 0 : CoordinateMaps::Distribution radial_distribution_envelope_ =
616 : CoordinateMaps::Distribution::Projective;
617 0 : std::vector<double> radial_partitioning_outer_shell_;
618 0 : std::vector<CoordinateMaps::Distribution> radial_distribution_outer_shell_ = {
619 : CoordinateMaps::Distribution::Linear};
620 0 : double translation_{};
621 0 : double length_inner_cube_{};
622 0 : double length_outer_cube_{};
623 0 : size_t number_of_outer_shells_{};
624 0 : size_t number_of_blocks_{};
625 0 : size_t first_outer_shell_block_{};
626 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
627 0 : outer_boundary_condition_;
628 0 : std::vector<std::string> block_names_;
629 : std::unordered_map<std::string, std::unordered_set<std::string>>
630 0 : block_groups_;
631 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
632 0 : grid_anchors_;
633 0 : double offset_x_coord_a_{};
634 0 : double offset_x_coord_b_{};
635 :
636 : // Variables to handle std::variant on Object A and B
637 0 : double x_coord_a_{};
638 0 : double x_coord_b_{};
639 0 : bool is_excised_a_ = false;
640 0 : bool is_excised_b_ = false;
641 0 : bool use_single_block_a_ = false;
642 0 : bool use_single_block_b_ = false;
643 0 : std::optional<bco::TimeDependentMapOptions<false>> time_dependent_options_;
644 0 : double opening_angle_ = std::numeric_limits<double>::signaling_NaN();
645 0 : bool spherical_harmonics_in_wavezone_ = false;
646 0 : bool use_worldtube_ = false;
647 : };
648 :
649 1 : namespace bco {
650 : /*!
651 : * \brief Create a set of centers of objects for the binary domains.
652 : *
653 : * \details Will add the following centers to the set:
654 : *
655 : * - Center: The origin
656 : * - CenterA: Center of object A
657 : * - CenterB: Center of object B
658 : *
659 : * \return Object required by the DomainCreator%s
660 : */
661 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
662 1 : create_grid_anchors(const std::array<double, 3>& center_a,
663 : const std::array<double, 3>& center_b);
664 :
665 : /*!
666 : * \brief Validate `InitialRefinement` map entries.
667 : *
668 : * \details Any spherical-harmonic block must use `size_t` (radial only).
669 : * `array<3>` is rejected on spherical-harmonic blocks even if angular
670 : * components are zero. Non-spherical-harmonic entries must use `array<3>`.
671 : *
672 : * \param context options context
673 : * \param initial_refinement the initial refinement from options
674 : * \param spherical_harmonic_shell_names the names of spherical shell blocks
675 : * or groups that use spherical harmonics
676 : * \param cylinder_names the names of cylindrical blocks or groups
677 : */
678 1 : void validate_initial_refinement(
679 : const Options::Context& context,
680 : const BinaryCompactObject::InitialRefinement::type& initial_refinement,
681 : const std::unordered_set<std::string>& spherical_harmonic_shell_names,
682 : const std::unordered_set<std::string>& cylinder_names = {});
683 :
684 : /*!
685 : * \brief Validate `InitialGridPoints` map entries.
686 : *
687 : * \details Any spherical-harmonic shell block must use
688 : * `array<2>{radial, L_max}`. Non-spherical-harmonic blocks must use `array<3>`.
689 : *
690 : * \param context options context
691 : * \param initial_number_of_grid_points the initial grid points from options
692 : * \param spherical_harmonic_shell_names the names of spherical shell blocks
693 : * or groups that use spherical harmonics
694 : * \param filled_cylinder_names the names of filled cylindrical blocks or groups
695 : * containing them
696 : */
697 1 : void validate_initial_grid_points(
698 : const Options::Context& context,
699 : const BinaryCompactObject::InitialGridPoints::type&
700 : initial_number_of_grid_points,
701 : const std::unordered_set<std::string>& spherical_harmonic_shell_names,
702 : const std::unordered_set<std::string>& filled_cylinder_names = {});
703 :
704 : /*!
705 : * \brief Convert `size_t` radial h refinement entries for spherical harmonic
706 : * blocks to `{r, 0, 0}`
707 : *
708 : * \details All `array<3>` entries are unchanged.
709 : *
710 : * \param expand_over_blocks `ExpandOverBlocks` containing the block names and
711 : * block groups
712 : * \param initial_refinement the initial refinement from options
713 : * \param spherical_harmonic_shell_names the names of spherical shell blocks
714 : * or groups that use spherical harmonics
715 : * \param cylinder_names the names of cylindrical blocks or groups
716 : *
717 : * \return converted refinement
718 : */
719 1 : std::vector<std::array<size_t, 3>> set_initial_refinement(
720 : const ExpandOverBlocks<std::array<size_t, 3>>& expand_over_blocks,
721 : const BinaryCompactObject::InitialRefinement::type& initial_refinement,
722 : const std::unordered_set<std::string>& spherical_harmonic_shell_names,
723 : const std::unordered_set<std::string>& cylinder_names = {});
724 :
725 : /*!
726 : * \brief Convert `array<2>{r, l_max}` entries for spherical harmonic blocks to
727 : * `{r, l_max, m_max}` = `{r, l_max, l_max}`
728 : *
729 : * \details All array<3> entries are used unchanged. We store the (l_max, m_max)
730 : * of the shell directly because ell is clear and unambiguous, whereas the
731 : * number of collocation points implied by ell depends on the spectral
732 : * implementation. The conversion to the number of collocation points is applied
733 : * in `initial_extents()`.
734 : *
735 : * \param expand_over_blocks `ExpandOverBlocks` containing the block names and
736 : * block groups
737 : * \param initial_grid_points the initial grid points from options
738 : * \param spherical_harmonic_shell_names the names of spherical shell blocks
739 : * or groups that use spherical harmonics
740 : * \param filled_cylinder_names the names of filled cylindrical blocks or groups
741 : * containing them
742 : *
743 : * \return converted grid points
744 : */
745 1 : std::vector<std::array<size_t, 3>> set_initial_grid_points(
746 : const ExpandOverBlocks<std::array<size_t, 3>>& expand_over_blocks,
747 : const BinaryCompactObject::InitialGridPoints::type& initial_grid_points,
748 : const std::unordered_set<std::string>& spherical_harmonic_shell_names,
749 : const std::unordered_set<std::string>& filled_cylinder_names = {});
750 : } // namespace bco
751 : } // namespace domain::creators
|
