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 Identity = CoordinateMaps::Identity<1>;
158 0 : using Identity2D = CoordinateMaps::Identity<2>;
159 0 : using Interval = CoordinateMaps::Interval;
160 0 : using RadialInterval3D =
161 : CoordinateMaps::ProductOf3Maps<Identity, Identity, Interval>;
162 : // The Translation type is no longer needed, but it is kept here for backwards
163 : // compatibility with old domains.
164 0 : using Translation = CoordinateMaps::ProductOf2Maps<Affine, Identity2D>;
165 0 : using Equiangular = CoordinateMaps::Equiangular;
166 0 : using Equiangular3D =
167 : CoordinateMaps::ProductOf3Maps<Equiangular, Equiangular, Equiangular>;
168 :
169 : public:
170 0 : using maps_list = tmpl::flatten<tmpl::list<
171 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Affine3D>,
172 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
173 : Equiangular3D>,
174 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Affine3D,
175 : Translation>,
176 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
177 : CoordinateMaps::DiscreteRotation<3>, Affine3D>,
178 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
179 : Equiangular3D>,
180 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Equiangular3D,
181 : Translation>,
182 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
183 : CoordinateMaps::Frustum>,
184 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
185 : CoordinateMaps::Wedge<3>>,
186 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
187 : CoordinateMaps::Wedge<3>, Translation>,
188 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Affine3D,
189 : Affine3D>,
190 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Equiangular3D,
191 : Affine3D>,
192 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
193 : CoordinateMaps::Wedge<3>, Affine3D>,
194 : domain::CoordinateMap<
195 : Frame::BlockLogical, Frame::Inertial,
196 : domain::CoordinateMaps::ProductOf2Maps<
197 : domain::CoordinateMaps::Interval,
198 : domain::CoordinateMaps::Identity<2>>,
199 : domain::CoordinateMaps::SphericalToCartesianPfaffian>,
200 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
201 : RadialInterval3D, CoordinateMaps::Wedge<3>,
202 : Affine3D>,
203 : bco::TimeDependentMapOptions<false>::maps_list>>;
204 :
205 : /// Options for an excision region in the domain
206 1 : struct Excision {
207 0 : static constexpr Options::String help = {
208 : "Excise the interior of the object, leaving a spherical hole in its "
209 : "absence."};
210 : template <typename BoundaryConditionsBase>
211 0 : struct BoundaryCondition {
212 0 : static std::string name() { return "ExciseWithBoundaryCondition"; }
213 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
214 0 : static constexpr Options::String help = {
215 : "The boundary condition to impose on the excision surface."};
216 : };
217 : template <typename Metavariables>
218 0 : using options = tmpl::list<BoundaryCondition<
219 : domain::BoundaryConditions::get_boundary_conditions_base<
220 : typename Metavariables::system>>>;
221 0 : Excision() = default;
222 : // NOLINTNEXTLINE(google-explicit-constructor)
223 0 : Excision(std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
224 : boundary_condition_in)
225 : : boundary_condition(std::move(boundary_condition_in)) {}
226 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
227 0 : boundary_condition;
228 : };
229 :
230 : /// Options for one of the two objects in the binary domain
231 1 : struct Object {
232 0 : static constexpr Options::String help = {
233 : "Options for an object in a binary domain."};
234 0 : struct InnerRadius {
235 0 : using type = double;
236 0 : static constexpr Options::String help = {
237 : "Inner coordinate radius of Layer 1."};
238 0 : static double lower_bound() { return 0.; }
239 : };
240 0 : struct OuterRadius {
241 0 : using type = double;
242 0 : static constexpr Options::String help = {
243 : "Outer coordinate radius of Layer 1"};
244 0 : static double lower_bound() { return 0.; }
245 : };
246 0 : struct XCoord {
247 0 : using type = double;
248 0 : static constexpr Options::String help = {"x-coordinate of center."};
249 : };
250 0 : struct Interior {
251 0 : using type = Options::Auto<Excision>;
252 0 : static constexpr Options::String help = {
253 : "Specify 'ExciseWithBoundaryCondition' and a boundary condition to "
254 : "excise Layer 0, leaving a spherical hole in its absence, or set to "
255 : "'Auto' to fill the interior."};
256 : };
257 0 : struct ExciseInterior {
258 0 : using type = bool;
259 0 : static constexpr Options::String help = {
260 : "Excise Layer 0, leaving a spherical hole in its absence."};
261 : };
262 0 : struct ShellLogMapStrength {
263 0 : using type = Options::Auto<double, Options::AutoLabel::None>;
264 0 : static constexpr Options::String help = {
265 : "Strength of the radial logarithmic map applied to the spherical "
266 : "shell enveloping the object. "
267 : "Set to 'None' for a linear radial grid. "
268 : "A value of 1 places the log-map origin (the 'singularity position') "
269 : "at radius 0. Larger values move the origin closer to the inner "
270 : "radius, concentrating grid points near the object. "
271 : "Enabling the log map requires the interior to be excised."};
272 : };
273 0 : struct CubeLogMapStrength {
274 0 : using type = Options::Auto<double, Options::AutoLabel::None>;
275 0 : static constexpr Options::String help = {
276 : "Strength of the radial logarithmic map applied to the cube "
277 : "surrounding the shell of each object. "
278 : "Set to 'None' for a linear radial grid. "
279 : "A value of 1 places the log-map origin (the 'singularity position') "
280 : "at radius 0. Larger values move the origin closer to the inner "
281 : "radius, concentrating grid points near the object."};
282 : };
283 : template <typename Metavariables>
284 0 : using options = tmpl::list<
285 : InnerRadius, OuterRadius, XCoord,
286 : tmpl::conditional_t<
287 : domain::BoundaryConditions::has_boundary_conditions_base_v<
288 : typename Metavariables::system>,
289 : Interior, ExciseInterior>,
290 : ShellLogMapStrength, CubeLogMapStrength>;
291 0 : Object() {} // NOLINT(modernize-use-equals-default)
292 0 : Object(double local_inner_radius, double local_outer_radius,
293 : double local_x_coord, std::optional<Excision> interior,
294 : std::optional<double> local_shell_log_map_strength,
295 : std::optional<double> local_cube_log_map_strength = std::nullopt)
296 : : inner_radius(local_inner_radius),
297 : outer_radius(local_outer_radius),
298 : x_coord(local_x_coord),
299 : inner_boundary_condition(
300 : interior.has_value()
301 : ? std::make_optional(std::move(interior->boundary_condition))
302 : : std::nullopt),
303 : shell_log_map_strength(local_shell_log_map_strength),
304 : cube_log_map_strength(local_cube_log_map_strength) {}
305 0 : Object(double local_inner_radius, double local_outer_radius,
306 : double local_x_coord, bool local_excise_interior,
307 : std::optional<double> local_shell_log_map_strength,
308 : std::optional<double> local_cube_log_map_strength = std::nullopt)
309 : : inner_radius(local_inner_radius),
310 : outer_radius(local_outer_radius),
311 : x_coord(local_x_coord),
312 : inner_boundary_condition(
313 : local_excise_interior
314 : ? std::optional<std::unique_ptr<
315 : domain::BoundaryConditions::BoundaryCondition>>{nullptr}
316 : : std::nullopt),
317 : shell_log_map_strength(local_shell_log_map_strength),
318 : cube_log_map_strength(local_cube_log_map_strength) {}
319 :
320 : /// Whether or not the object should be excised from the domain, leaving a
321 : /// spherical hole. When this is true, `inner_boundary_condition` is
322 : /// guaranteed to hold a value (though it might be a `nullptr` if we are not
323 : /// working with boundary conditions).
324 1 : bool is_excised() const;
325 :
326 0 : double inner_radius{};
327 0 : double outer_radius{};
328 0 : double x_coord{};
329 : std::optional<
330 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>
331 0 : inner_boundary_condition;
332 0 : std::optional<double> shell_log_map_strength{};
333 0 : std::optional<double> cube_log_map_strength{};
334 : };
335 :
336 : // Simpler version of an object: a single cube centered on (xCoord,0,0)
337 0 : struct CartesianCubeAtXCoord {
338 0 : static constexpr Options::String help = {
339 : "Options to set a single cube at a location on the x-axis"};
340 0 : struct XCoord {
341 0 : static std::string name() { return "CartesianCubeAtXCoord"; }
342 0 : using type = double;
343 0 : static constexpr Options::String help = {"x-coordinate of center."};
344 : };
345 0 : using options = tmpl::list<XCoord>;
346 0 : CartesianCubeAtXCoord() = default;
347 : // NOLINTNEXTLINE(google-explicit-constructor)
348 0 : CartesianCubeAtXCoord(const double x_coord_in) : x_coord(x_coord_in) {}
349 0 : bool is_excised() const { return false; }
350 0 : double x_coord;
351 : };
352 :
353 0 : struct ObjectA {
354 0 : using type = std::variant<Object, CartesianCubeAtXCoord>;
355 0 : static constexpr Options::String help = {
356 : "Options for the object to the right of the origin (along the positive "
357 : "x-axis)."};
358 : };
359 :
360 0 : struct ObjectB {
361 0 : using type = std::variant<Object, CartesianCubeAtXCoord>;
362 0 : static constexpr Options::String help = {
363 : "Options for the object to the left of the origin (along the negative "
364 : "x-axis)."};
365 : };
366 :
367 0 : struct CenterOfMassOffset {
368 0 : using type = std::array<double, 2>;
369 0 : static constexpr Options::String help = {
370 : "Offset in the y and z axes applied to both object A and B in order to "
371 : "control the center of mass. This moves the location of the two objects"
372 : " in the grid frame but keeps the Envelope and OuterShell centered on "
373 : "the origin in the grid frame."};
374 : };
375 :
376 0 : struct Envelope {
377 0 : static constexpr Options::String help = {
378 : "Options for the sphere enveloping the two objects."};
379 : };
380 :
381 0 : struct EnvelopeRadius {
382 0 : using group = Envelope;
383 0 : static std::string name() { return "Radius"; }
384 0 : using type = double;
385 0 : static constexpr Options::String help = {
386 : "Radius of the sphere enveloping the two objects."};
387 : };
388 :
389 0 : struct OuterShell {
390 0 : static constexpr Options::String help = {
391 : "Options for the outer spherical shell."};
392 : };
393 :
394 0 : struct OuterRadius {
395 0 : using group = OuterShell;
396 0 : static std::string name() { return "Radius"; }
397 0 : using type = double;
398 0 : static constexpr Options::String help = {"Radius of the entire domain."};
399 : };
400 :
401 0 : struct RadialPartitioningOuterShell {
402 0 : static std::string name() { return "RadialPartitioning"; }
403 0 : using group = OuterShell;
404 0 : using type = std::vector<double>;
405 0 : static constexpr Options::String help = {
406 : "Radial coordinates of the boundaries splitting the outer spherical "
407 : "shells between the envelope radius and OuterRadius. They must be "
408 : "given in ascending order. This should be used if boundaries need to "
409 : "be set at specific radii. If the number but not the specific "
410 : "locations of the boundaries are important, use InitialRefinement "
411 : "instead."};
412 : };
413 :
414 0 : struct RadialDistributionOuterShell {
415 0 : static std::string name() { return "RadialDistribution"; }
416 0 : using group = OuterShell;
417 0 : using type =
418 : std::variant<domain::CoordinateMaps::Distribution,
419 : std::vector<domain::CoordinateMaps::Distribution>>;
420 0 : static constexpr Options::String help = {
421 : "Select the radial distribution of grid points in each outer spherical "
422 : "shell. There must be N+1 radial distributions specified for N radial "
423 : "partitions. You can also specify just a single radial distribution "
424 : "(not in a vector) which will use the same distribution for all "
425 : "partitions."};
426 : };
427 :
428 0 : struct OpeningAngle {
429 0 : using group = OuterShell;
430 0 : static std::string name() { return "OpeningAngle"; }
431 0 : using type = double;
432 0 : static constexpr Options::String help = {
433 : "The combined opening angle of the two half wedges of the outer shell"
434 : " in degrees. A value of 120.0 partitions the x-y and x-z slices of the"
435 : " outer shell into six Blocks of equal angular size."};
436 : };
437 :
438 0 : struct SphericalHarmonicsInWavezone {
439 0 : using group = OuterShell;
440 0 : static std::string name() { return "UseSphericalHarmonics"; }
441 0 : using type = bool;
442 0 : static bool suggested_value() { return false; }
443 0 : static constexpr Options::String help = {
444 : "Use a spherical-harmonic basis for the outer wavezone shell(s) "
445 : "instead of the default 10-wedge Cartesian basis. When enabled, "
446 : "InitialGridPoints for wavezone blocks must be specified as "
447 : "array<size_t, 2> = {radial_points, L_max}, and InitialRefinement "
448 : "as a scalar or single-element value (angular refinement is fixed at "
449 : "0). "};
450 : };
451 :
452 0 : struct UseWorldtube {
453 0 : using type = bool;
454 0 : static constexpr Options::String help = {
455 : "Whether to set up functions of time appropriate for a worldtube run. "
456 : "When true, some functions of time will be IntegratedFunctionOfTime "
457 : "used to control the orbit of the worldtube. Used by the curved scalar "
458 : "wave worldtube executable."};
459 : };
460 :
461 0 : struct CubeScale {
462 0 : using type = double;
463 0 : static constexpr Options::String help = {
464 : "Specify the desired cube scale that must be greater than or equal to "
465 : "1.0. The initial separation is multiplied by this cube scale to "
466 : "produce larger cubes around each object which is desirable when "
467 : "closer to merger."};
468 0 : static double lower_bound() { return 1.0; }
469 : };
470 :
471 0 : struct InitialRefinement {
472 0 : using type = std::variant<
473 : size_t, std::array<size_t, 3>, std::vector<std::array<size_t, 3>>,
474 : std::unordered_map<std::string,
475 : std::variant<std::array<size_t, 3>, size_t>>>;
476 0 : static constexpr Options::String help = {
477 : "Initial refinement level in each block of the domain. See main help "
478 : "text for details."};
479 : };
480 :
481 0 : struct InitialGridPoints {
482 0 : using type = std::variant<
483 : size_t, std::array<size_t, 3>, std::vector<std::array<size_t, 3>>,
484 : std::unordered_map<std::string, std::array<size_t, 3>>,
485 : std::unordered_map<std::string, std::variant<std::array<size_t, 3>,
486 : std::array<size_t, 2>>>>;
487 0 : static constexpr Options::String help = {
488 : "Initial number of grid points in the elements of each block of the "
489 : "domain. See main help text for details."};
490 : };
491 :
492 0 : struct UseEquiangularMap {
493 0 : using type = bool;
494 0 : static constexpr Options::String help = {
495 : "Distribute grid points equiangularly."};
496 0 : static bool suggested_value() { return true; }
497 : };
498 :
499 0 : struct RadialDistributionEnvelope {
500 0 : using group = Envelope;
501 0 : static std::string name() { return "RadialDistribution"; }
502 0 : using type = CoordinateMaps::Distribution;
503 0 : static constexpr Options::String help = {
504 : "The distribution of radial grid points in the envelope, the layer "
505 : "made of ten bulged Frustums."};
506 : };
507 :
508 : template <typename BoundaryConditionsBase>
509 0 : struct OuterBoundaryCondition {
510 0 : using group = OuterShell;
511 0 : static std::string name() { return "BoundaryCondition"; }
512 0 : static constexpr Options::String help =
513 : "Options for the outer boundary conditions.";
514 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
515 : };
516 :
517 : // This is for optional time dependent maps
518 0 : struct TimeDependentMaps {
519 0 : using type = Options::Auto<bco::TimeDependentMapOptions<false>,
520 : Options::AutoLabel::None>;
521 0 : static constexpr Options::String help =
522 : bco::TimeDependentMapOptions<false>::help;
523 : };
524 :
525 : template <typename Metavariables>
526 0 : using options = tmpl::append<
527 : tmpl::list<ObjectA, ObjectB, CenterOfMassOffset, EnvelopeRadius,
528 : OuterRadius, CubeScale, InitialRefinement, InitialGridPoints,
529 : UseEquiangularMap, RadialDistributionEnvelope,
530 : RadialPartitioningOuterShell, RadialDistributionOuterShell,
531 : OpeningAngle, SphericalHarmonicsInWavezone, UseWorldtube,
532 : TimeDependentMaps>,
533 : tmpl::conditional_t<
534 : domain::BoundaryConditions::has_boundary_conditions_base_v<
535 : typename Metavariables::system>,
536 : tmpl::list<OuterBoundaryCondition<
537 : domain::BoundaryConditions::get_boundary_conditions_base<
538 : typename Metavariables::system>>>,
539 : tmpl::list<>>>;
540 :
541 0 : static constexpr Options::String help{
542 : "A general domain for two compact objects. Each object is represented by "
543 : "a cube along the x-axis. Object A is located on the right and Object B "
544 : "is located on the left. Their locations should be chosen such that "
545 : "their center of mass is located at the origin."
546 : "The interior of each object can have a spherical excision to "
547 : "represent a black hole."
548 : "\n"
549 : "The two objects are enveloped by a sphere centered at the origin, "
550 : "and by an outer shell that can transition to large outer radii."
551 : "\n"
552 : "Both the InitialRefinement and the InitialGridPoints can be one of "
553 : "the following:\n"
554 : " - A single number: Uniform refinement in all blocks and "
555 : "dimensions\n"
556 : " - Three numbers: Refinement in [polar, azimuthal, radial] direction "
557 : "in all blocks\n"
558 : " - A map from block names or groups to three numbers: Per-block "
559 : "refinement in [polar, azimuthal, radial] direction\n"
560 : " - A list, with [polar, azimuthal, radial] refinement for each block\n"
561 : "\n"
562 : "The domain can rotate around the "
563 : "z-axis and expand/compress radially. The two objects can each have a "
564 : "spherical distortion (shape map)."};
565 :
566 0 : BinaryCompactObject(
567 : typename ObjectA::type object_A, typename ObjectB::type object_B,
568 : std::array<double, 2> center_of_mass_offset, double envelope_radius,
569 : double outer_radius, double cube_scale,
570 : const typename InitialRefinement::type& initial_refinement,
571 : const typename InitialGridPoints::type& initial_number_of_grid_points,
572 : bool use_equiangular_map = true,
573 : CoordinateMaps::Distribution radial_distribution_envelope =
574 : CoordinateMaps::Distribution::Projective,
575 : const std::vector<double>& radial_partitioning_outer_shell = {},
576 : const typename RadialDistributionOuterShell::type&
577 : radial_distribution_outer_shell =
578 : CoordinateMaps::Distribution::Linear,
579 : double opening_angle_in_degrees = 90.0,
580 : bool spherical_harmonics_in_wavezone = false, bool use_worldtube = false,
581 : std::optional<bco::TimeDependentMapOptions<false>>
582 : time_dependent_options = std::nullopt,
583 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
584 : outer_boundary_condition = nullptr,
585 : const Options::Context& context = {});
586 :
587 0 : BinaryCompactObject() = default;
588 0 : BinaryCompactObject(const BinaryCompactObject&) = delete;
589 0 : BinaryCompactObject(BinaryCompactObject&&) = default;
590 0 : BinaryCompactObject& operator=(const BinaryCompactObject&) = delete;
591 0 : BinaryCompactObject& operator=(BinaryCompactObject&&) = default;
592 0 : ~BinaryCompactObject() override = default;
593 :
594 0 : Domain<3> create_domain() const override;
595 :
596 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
597 1 : grid_anchors() const override {
598 : return grid_anchors_;
599 : }
600 :
601 : std::vector<DirectionMap<
602 : 3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
603 1 : external_boundary_conditions() const override;
604 :
605 1 : std::vector<std::array<size_t, 3>> initial_extents() const override;
606 :
607 1 : std::vector<std::array<size_t, 3>> initial_refinement_levels()
608 : const override {
609 : return initial_refinement_;
610 : }
611 :
612 1 : std::vector<std::string> block_names() const override { return block_names_; }
613 :
614 : std::unordered_map<std::string, std::unordered_set<std::string>>
615 1 : block_groups() const override {
616 : return block_groups_;
617 : }
618 :
619 1 : auto functions_of_time(const std::unordered_map<std::string, double>&
620 : initial_expiration_times = {}) const
621 : -> std::unordered_map<
622 : std::string,
623 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>> override;
624 :
625 : private:
626 0 : typename ObjectA::type object_A_{Object{}};
627 0 : typename ObjectB::type object_B_{Object{}};
628 0 : std::array<double, 2> center_of_mass_offset_{};
629 0 : double envelope_radius_ = std::numeric_limits<double>::signaling_NaN();
630 0 : double outer_radius_ = std::numeric_limits<double>::signaling_NaN();
631 0 : std::vector<std::array<size_t, 3>> initial_refinement_;
632 : // For most blocks this stores the number of grid points in each logical
633 : // direction. For spherical-harmonic shell blocksthe angular entries instead
634 : // store the spherical-harmonic degrees {n_radial, l_max, m_max}, because
635 : // ell is unambiguous whereas the implied number of collocation points
636 : // depends on the spectral implementation. The conversion of (l_max, m_max)
637 : // to the actual number of collocation points is applied in
638 : // `initial_extents()`.
639 0 : std::vector<std::array<size_t, 3>> initial_number_of_grid_points_;
640 0 : bool use_equiangular_map_ = true;
641 0 : CoordinateMaps::Distribution radial_distribution_envelope_ =
642 : CoordinateMaps::Distribution::Projective;
643 0 : std::vector<double> radial_partitioning_outer_shell_;
644 0 : std::vector<CoordinateMaps::Distribution> radial_distribution_outer_shell_ = {
645 : CoordinateMaps::Distribution::Linear};
646 0 : double translation_{};
647 0 : double length_inner_cube_{};
648 0 : double length_outer_cube_{};
649 0 : size_t number_of_outer_shells_{};
650 0 : size_t number_of_blocks_{};
651 0 : size_t first_outer_shell_block_{};
652 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
653 0 : outer_boundary_condition_;
654 0 : std::vector<std::string> block_names_;
655 : std::unordered_map<std::string, std::unordered_set<std::string>>
656 0 : block_groups_;
657 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
658 0 : grid_anchors_;
659 0 : double offset_x_coord_a_{};
660 0 : double offset_x_coord_b_{};
661 :
662 : // Variables to handle std::variant on Object A and B
663 0 : double x_coord_a_{};
664 0 : double x_coord_b_{};
665 0 : bool is_excised_a_ = false;
666 0 : bool is_excised_b_ = false;
667 0 : bool use_single_block_a_ = false;
668 0 : bool use_single_block_b_ = false;
669 0 : std::optional<bco::TimeDependentMapOptions<false>> time_dependent_options_;
670 0 : double opening_angle_ = std::numeric_limits<double>::signaling_NaN();
671 0 : bool spherical_harmonics_in_wavezone_ = false;
672 0 : bool use_worldtube_ = false;
673 : };
674 :
675 1 : namespace bco {
676 : /// Superset of the InitialRefinement option types for the variants used by
677 : /// `domain::creators::BinaryCompactObject` and
678 : /// `domain::creators::CylindricalBinaryCompactObject`, which is useful for
679 : /// reusing functionality, e.g. `bco::validate_initial_refinement`
680 1 : using InitialRefinement = std::variant<
681 : size_t, std::array<size_t, 3>, std::vector<std::array<size_t, 3>>,
682 : std::unordered_map<std::string,
683 : std::variant<std::array<size_t, 3>, size_t>>,
684 : std::unordered_map<std::string,
685 : std::variant<std::array<size_t, 2>, size_t>>>;
686 :
687 : /*!
688 : * \brief Create a set of centers of objects for the binary domains.
689 : *
690 : * \details Will add the following centers to the set:
691 : *
692 : * - Center: The origin
693 : * - CenterA: Center of object A
694 : * - CenterB: Center of object B
695 : *
696 : * \return Object required by the DomainCreator%s
697 : */
698 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
699 1 : create_grid_anchors(const std::array<double, 3>& center_a,
700 : const std::array<double, 3>& center_b);
701 :
702 : /*!
703 : * \brief Validate `InitialRefinement` map entries.
704 : *
705 : * \details Spherical-harmonic blocks must use `size_t` (radial only),
706 : * cylindrical blocks must use `array<2>{radial, z}`, and all other blocks must
707 : * use `array<3>`.
708 : *
709 : * \param context options context
710 : * \param initial_refinement the initial refinement from options
711 : * \param spherical_harmonic_shell_names the names of spherical shell blocks
712 : * or groups that use spherical harmonics
713 : * \param cylinder_names the names of cylindrical blocks or groups
714 : */
715 1 : void validate_initial_refinement(
716 : const Options::Context& context,
717 : const InitialRefinement& initial_refinement,
718 : const std::unordered_set<std::string>& spherical_harmonic_shell_names,
719 : const std::unordered_set<std::string>& cylinder_names = {});
720 :
721 : /*!
722 : * \brief Validate `InitialGridPoints` map entries.
723 : *
724 : * \details Any spherical-harmonic shell block must use
725 : * `array<2>{radial, L_max}`. Non-spherical-harmonic blocks must use `array<3>`.
726 : *
727 : * \param context options context
728 : * \param initial_number_of_grid_points the initial grid points from options
729 : * \param spherical_harmonic_shell_names the names of spherical shell blocks
730 : * or groups that use spherical harmonics
731 : * \param hollow_cylinder_names the names of hollow cylindrical blocks
732 : * \param filled_cylinder_names the names of filled cylindrical blocks or groups
733 : * containing them
734 : */
735 1 : void validate_initial_grid_points(
736 : const Options::Context& context,
737 : const BinaryCompactObject::InitialGridPoints::type&
738 : initial_number_of_grid_points,
739 : const std::unordered_set<std::string>& spherical_harmonic_shell_names,
740 : const std::unordered_set<std::string>& hollow_cylinder_names = {},
741 : const std::unordered_set<std::string>& filled_cylinder_names = {});
742 :
743 : /*!
744 : * \brief Convert `size_t` radial h refinement entries for spherical harmonic
745 : * blocks to `{r, 0, 0}` and `array<2>`[r, z] radial and z h refinement entries
746 : * for cylindrical blocks to `{r, 0, z}`
747 : *
748 : * \details All `array<3>` entries are unchanged.
749 : *
750 : * \param expand_over_blocks `ExpandOverBlocks` containing the block names and
751 : * block groups
752 : * \param initial_refinement the initial refinement from options
753 : * \return converted refinement
754 : */
755 1 : std::vector<std::array<size_t, 3>> set_initial_refinement(
756 : const ExpandOverBlocks<std::array<size_t, 3>>& expand_over_blocks,
757 : const InitialRefinement& initial_refinement);
758 :
759 : /*!
760 : * \brief Convert `array<2>{r, l_max}` entries for spherical harmonic blocks to
761 : * `{r, l_max, m_max}` = `{r, l_max, l_max}`
762 : *
763 : * \details All array<3> entries are used unchanged. We store the (l_max, m_max)
764 : * of the shell directly because ell is clear and unambiguous, whereas the
765 : * number of collocation points implied by ell depends on the spectral
766 : * implementation. The conversion to the number of collocation points is applied
767 : * in `initial_extents()`.
768 : *
769 : * \param expand_over_blocks `ExpandOverBlocks` containing the block names and
770 : * block groups
771 : * \param initial_grid_points the initial grid points from options
772 : * \param spherical_harmonic_shell_names the names of spherical shell blocks
773 : * or groups that use spherical harmonics
774 : * \param filled_cylinder_names the names of filled cylindrical blocks or groups
775 : * containing them
776 : *
777 : * \return converted grid points
778 : */
779 1 : std::vector<std::array<size_t, 3>> set_initial_grid_points(
780 : const ExpandOverBlocks<std::array<size_t, 3>>& expand_over_blocks,
781 : const BinaryCompactObject::InitialGridPoints::type& initial_grid_points,
782 : const std::unordered_set<std::string>& spherical_harmonic_shell_names,
783 : const std::unordered_set<std::string>& filled_cylinder_names = {});
784 : } // namespace bco
785 : } // namespace domain::creators
|
