Skip to content

Commit

Permalink
Add BCO Offset Option
Browse files Browse the repository at this point in the history 
Co-authored-by: Alexandra Macedo <alexandra.l.macedo@gmail.com>
Co-authored-by: Marceline Bonilla <marceline.bonilla@black-holes.org>
  • Loading branch information
@AlexCarpenter46 @macedo22 @mar-celine
3 people committed Sep 6, 2024
1 parent 27756af commit 66fb37d
Show file tree
Hide file tree
Showing 26 changed files with 852 additions and 451 deletions.
Binary file modified docs/Images/Domain/DomainCreators/binary_compact_object_domain.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
153 changes: 111 additions & 42 deletions src/Domain/Creators/BinaryCompactObject.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@
#include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp"
#include "Domain/Structure/BlockNeighbor.hpp"
#include "Options/ParseError.hpp"
#include "Utilities/EqualWithinRoundoff.hpp"
#include "Utilities/MakeArray.hpp"

namespace Frame {
Expand Down Expand Up @@ -73,7 +74,7 @@ template <bool UseWorldtube>
BinaryCompactObject<UseWorldtube>::BinaryCompactObject(
typename ObjectA::type object_A, typename ObjectB::type object_B,
std::array<double, 2> center_of_mass_offset, const double envelope_radius,
const double outer_radius,
const double outer_radius, const double cube_scale,
const typename InitialRefinement::type& initial_refinement,
const typename InitialGridPoints::type& initial_number_of_grid_points,
const bool use_equiangular_map,
Expand Down Expand Up @@ -113,10 +114,31 @@ BinaryCompactObject<UseWorldtube>::BinaryCompactObject(
// Determination of parameters for domain construction:
const double tan_half_opening_angle = tan(0.5 * opening_angle_);
translation_ = 0.5 * (x_coord_a_ + x_coord_b_);
length_inner_cube_ = abs(x_coord_a_ - x_coord_b_);
length_inner_cube_ = cube_scale * (x_coord_a_ - x_coord_b_);
if (length_inner_cube_ < (x_coord_a_ - x_coord_b_)) {
PARSE_ERROR(
context,
"The cube length should be greater than or equal to the initial "
"separation between the two objects.");
}
length_outer_cube_ =
2.0 * envelope_radius_ / sqrt(2.0 + square(tan_half_opening_angle));

// We chose to handle potential roundoff differences here by using equal
// within roundoff instead of exact equality because the wedge map expects
// exact equality when checking for a zero offset.
if (equal_within_roundoff(length_inner_cube_, (x_coord_a_ - x_coord_b_),
std::numeric_limits<double>::epsilon() * 100.0,
length_inner_cube_)) {
offset_x_coord_a_ = 0.0;
offset_x_coord_b_ = 0.0;
} else {
offset_x_coord_a_ =
x_coord_a_ - (x_coord_a_ + x_coord_b_ + length_inner_cube_) * 0.5;
offset_x_coord_b_ =
x_coord_b_ - (x_coord_a_ + x_coord_b_ - length_inner_cube_) * 0.5;
}

// Calculate number of blocks
// Object cubes and shells have 6 blocks each, for a total for 24 blocks.
// The envelope and outer shell have another 10 blocks each.
Expand Down Expand Up @@ -202,6 +224,16 @@ BinaryCompactObject<UseWorldtube>::BinaryCompactObject(
"or neither.");
}
}
const bool filled_excision_a = not(use_single_block_a_ or is_excised_a_);
const bool filled_excision_b = not(use_single_block_b_ or is_excised_b_);
if ((filled_excision_a or filled_excision_b) and
not equal_within_roundoff(offset_x_coord_a_, 0.0)) {
PARSE_ERROR(context,
"Setting CubeScale > 1.0 is not supported for domains with "
"ExciseInterior = False. Consider using "
"CartesianCubeAtX for the Object without an excised interior.");
}

if (envelope_radius_ >= outer_radius_) {
PARSE_ERROR(context,
"The outer radius must be larger than the envelope radius.");
Expand Down Expand Up @@ -329,12 +361,25 @@ BinaryCompactObject<UseWorldtube>::BinaryCompactObject(
}

if (time_dependent_options_.has_value()) {
// The reason we don't just always use half the inner cube length is to
// avoid potential roundoff issues if there is no offset
const std::optional<std::array<double, 3>> cube_A_center =
length_inner_cube_ == x_coord_a_ - x_coord_b_
? std::optional<std::array<double, 3>>{}
: std::array{translation_ + 0.5 * length_inner_cube_,
center_of_mass_offset_[0], center_of_mass_offset_[1]};
const std::optional<std::array<double, 3>> cube_B_center =
length_inner_cube_ == x_coord_a_ - x_coord_b_
? std::optional<std::array<double, 3>>{}
: std::array{translation_ - 0.5 * length_inner_cube_,
center_of_mass_offset_[0], center_of_mass_offset_[1]};
time_dependent_options_->build_maps(
std::array{std::array{x_coord_a_, center_of_mass_offset_[0],
center_of_mass_offset_[1]},
std::array{x_coord_b_, center_of_mass_offset_[0],
center_of_mass_offset_[1]}},
radii_A, radii_B, envelope_radius_, outer_radius_);
cube_A_center, cube_B_center, radii_A, radii_B, envelope_radius_,
outer_radius_);
}
}

Expand Down Expand Up @@ -364,13 +409,15 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {

// ObjectA/B is on the right/left, respectively.
const Affine3D translation_A{
Affine{-1.0, 1.0, -1.0 + x_coord_a_, 1.0 + x_coord_a_},
Affine{-1.0, 1.0, -1.0 + x_coord_a_ - offset_x_coord_a_,
1.0 + x_coord_a_ - offset_x_coord_a_},
Affine{-1.0, 1.0, -1.0 + center_of_mass_offset_[0],
1.0 + center_of_mass_offset_[0]},
Affine{-1.0, 1.0, -1.0 + center_of_mass_offset_[1],
1.0 + center_of_mass_offset_[1]}};
const Affine3D translation_B{
Affine{-1.0, 1.0, -1.0 + x_coord_b_, 1.0 + x_coord_b_},
Affine{-1.0, 1.0, -1.0 + x_coord_b_ - offset_x_coord_b_,
1.0 + x_coord_b_ - offset_x_coord_b_},
Affine{-1.0, 1.0, -1.0 + center_of_mass_offset_[0],
1.0 + center_of_mass_offset_[0]},
Affine{-1.0, 1.0, -1.0 + center_of_mass_offset_[1],
Expand All @@ -380,8 +427,9 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
if (use_single_block_a_) {
maps.emplace_back(
make_coordinate_map_base<Frame::BlockLogical, Frame::Inertial>(Affine3D{
Affine(-1.0, 1.0, -0.5 * length_inner_cube_ + x_coord_a_,
0.5 * length_inner_cube_ + x_coord_a_),
Affine(-1.0, 1.0,
-0.5 * length_inner_cube_ + x_coord_a_ - offset_x_coord_a_,
0.5 * length_inner_cube_ + x_coord_a_ - offset_x_coord_a_),
Affine(-1.0, 1.0,
-0.5 * length_inner_cube_ + center_of_mass_offset_[0],
0.5 * length_inner_cube_ + center_of_mass_offset_[0]),
Expand All @@ -394,21 +442,26 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
// Each object is surrounded by 6 inner wedges that make a sphere, and
// another 6 outer wedges that transition to a cube.
const auto& object_a = std::get<Object>(object_A_);

const auto& offset_a_optional =
offset_x_coord_a_ == 0
? std::nullopt
: std::make_optional(std::make_pair(
length_inner_cube_ * 0.5,
std::array<double, 3>{{offset_x_coord_a_, 0.0, 0.0}}));
Maps maps_center_A =
domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(object_a.inner_radius,
object_a.outer_radius, inner_sphericity_A,
1.0, use_equiangular_map_, false, {},
object_A_radial_distribution),
sph_wedge_coordinate_maps(
object_a.inner_radius, object_a.outer_radius,
inner_sphericity_A, 1.0, use_equiangular_map_,
offset_a_optional, false, {}, object_A_radial_distribution),
translation_A);
Maps maps_cube_A =
domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(object_a.outer_radius,
sqrt(3.0) * 0.5 * length_inner_cube_, 1.0,
0.0, use_equiangular_map_),
sph_wedge_coordinate_maps(
object_a.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_,
1.0, 0.0, use_equiangular_map_, offset_a_optional),
translation_A);
std::move(maps_center_A.begin(), maps_center_A.end(),
std::back_inserter(maps));
Expand All @@ -417,8 +470,9 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
if (use_single_block_b_) {
maps.emplace_back(
make_coordinate_map_base<Frame::BlockLogical, Frame::Inertial>(Affine3D{
Affine(-1.0, 1.0, -0.5 * length_inner_cube_ + x_coord_b_,
0.5 * length_inner_cube_ + x_coord_b_),
Affine(-1.0, 1.0,
-0.5 * length_inner_cube_ + x_coord_b_ - offset_x_coord_b_,
0.5 * length_inner_cube_ + x_coord_b_ - offset_x_coord_b_),
Affine(-1.0, 1.0,
-0.5 * length_inner_cube_ + center_of_mass_offset_[0],
0.5 * length_inner_cube_ + center_of_mass_offset_[0]),
Expand All @@ -431,20 +485,26 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
// Each object is surrounded by 6 inner wedges that make a sphere, and
// another 6 outer wedges that transition to a cube.
const auto& object_b = std::get<Object>(object_B_);
const auto& offset_b_optional =
offset_x_coord_b_ == 0
? std::nullopt
: std::make_optional(std::make_pair(
length_inner_cube_ * 0.5,
std::array<double, 3>{{offset_x_coord_b_, 0.0, 0.0}}));
Maps maps_center_B =
domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(object_b.inner_radius,
object_b.outer_radius, inner_sphericity_B,
1.0, use_equiangular_map_, false, {},
object_B_radial_distribution),
sph_wedge_coordinate_maps(
object_b.inner_radius, object_b.outer_radius,
inner_sphericity_B, 1.0, use_equiangular_map_,
offset_b_optional, false, {}, object_B_radial_distribution),
translation_B);
Maps maps_cube_B =
domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(object_b.outer_radius,
sqrt(3.0) * 0.5 * length_inner_cube_, 1.0,
0.0, use_equiangular_map_),
sph_wedge_coordinate_maps(
object_b.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_,
1.0, 0.0, use_equiangular_map_, offset_b_optional),
translation_B);
std::move(maps_center_B.begin(), maps_center_B.end(),
std::back_inserter(maps));
Expand Down Expand Up @@ -472,10 +532,13 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
std::back_inserter(maps));

// --- Outer spherical shell (10 blocks) ---
Maps maps_outer_shell = domain::make_vector_coordinate_map_base<
Frame::BlockLogical, Frame::Inertial, 3>(sph_wedge_coordinate_maps(
envelope_radius_, outer_radius_, 1.0, 1.0, use_equiangular_map_, true, {},
{radial_distribution_outer_shell_}, ShellWedges::All, opening_angle_));
Maps maps_outer_shell =
domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(envelope_radius_, outer_radius_, 1.0, 1.0,
use_equiangular_map_, std::nullopt, true,
{}, {radial_distribution_outer_shell_},
ShellWedges::All, opening_angle_));
std::move(maps_outer_shell.begin(), maps_outer_shell.end(),
std::back_inserter(maps));

Expand All @@ -492,18 +555,21 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
if (use_equiangular_map_) {
maps.emplace_back(
make_coordinate_map_base<Frame::BlockLogical, Frame::Inertial>(
Equiangular3D{Equiangular(-1.0, 1.0, -1.0 * scaled_r_inner_A,
scaled_r_inner_A),
Equiangular(-1.0, 1.0, -1.0 * scaled_r_inner_A,
scaled_r_inner_A),
Equiangular(-1.0, 1.0, -1.0 * scaled_r_inner_A,
scaled_r_inner_A)},
Equiangular3D{
Equiangular(-1.0, 1.0,
-1.0 * scaled_r_inner_A + offset_x_coord_a_,
scaled_r_inner_A + offset_x_coord_a_),
Equiangular(-1.0, 1.0, -1.0 * scaled_r_inner_A,
scaled_r_inner_A),
Equiangular(-1.0, 1.0, -1.0 * scaled_r_inner_A,
scaled_r_inner_A)},
translation_A));
} else {
maps.emplace_back(make_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial>(
Affine3D{
Affine(-1.0, 1.0, -1.0 * scaled_r_inner_A, scaled_r_inner_A),
Affine(-1.0, 1.0, -1.0 * scaled_r_inner_A + offset_x_coord_a_,
scaled_r_inner_A + offset_x_coord_a_),
Affine(-1.0, 1.0, -1.0 * scaled_r_inner_A, scaled_r_inner_A),
Affine(-1.0, 1.0, -1.0 * scaled_r_inner_A, scaled_r_inner_A)},
translation_A));
Expand Down Expand Up @@ -535,18 +601,21 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
if (use_equiangular_map_) {
maps.emplace_back(
make_coordinate_map_base<Frame::BlockLogical, Frame::Inertial>(
Equiangular3D{Equiangular(-1.0, 1.0, -1.0 * scaled_r_inner_B,
scaled_r_inner_B),
Equiangular(-1.0, 1.0, -1.0 * scaled_r_inner_B,
scaled_r_inner_B),
Equiangular(-1.0, 1.0, -1.0 * scaled_r_inner_B,
scaled_r_inner_B)},
Equiangular3D{
Equiangular(-1.0, 1.0,
-1.0 * scaled_r_inner_B + offset_x_coord_b_,
scaled_r_inner_B + offset_x_coord_b_),
Equiangular(-1.0, 1.0, -1.0 * scaled_r_inner_B,
scaled_r_inner_B),
Equiangular(-1.0, 1.0, -1.0 * scaled_r_inner_B,
scaled_r_inner_B)},
translation_B));
} else {
maps.emplace_back(make_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial>(
Affine3D{
Affine(-1.0, 1.0, -1.0 * scaled_r_inner_B, scaled_r_inner_B),
Affine(-1.0, 1.0, -1.0 * scaled_r_inner_B + offset_x_coord_b_,
scaled_r_inner_B + offset_x_coord_b_),
Affine(-1.0, 1.0, -1.0 * scaled_r_inner_B, scaled_r_inner_B),
Affine(-1.0, 1.0, -1.0 * scaled_r_inner_B, scaled_r_inner_B)},
translation_B));
Expand Down
30 changes: 25 additions & 5 deletions src/Domain/Creators/BinaryCompactObject.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,9 +123,17 @@ create_grid_anchors(const std::array<double, 3>& center_a,
* hemispheres. The cutting plane always intersects the x-axis at the origin.
* - The x-coordinate locations of the two objects should be chosen such that
* the center of mass is located at x=0.
* - Alternatively, one can replace the inner shell and cube blocks of each
* object with a single cartesian cube. This is less efficient, but allows
* testing of methods only coded on cartesian grids.
* - The cubes are first constructed at the origin. Then, they are translated
* left/right by their Object's x-coordinate and offset depending on the cube
* length.
* - The CubeScale option describes how to scale the length of the cube
* surrounding object A/B. It must be greater than or equal to 1.0 with 1.0
* meaning the side length of the cube is the initial physical separation
* between the two objects. If CubeScale is greater than 1.0, the centers of
* the two objects will be offset relative to the centers of the cubes.
* - Alternatively, one can replace the inner shell and cube blocks of each
* object with a single cartesian cube. This is less efficient, but allows
* testing of methods only coded on cartesian grids.
*
* \par Time dependence:
* The following time-dependent maps are applied:
Expand Down Expand Up @@ -376,6 +384,16 @@ class BinaryCompactObject : public DomainCreator<3> {
" outer shell into six Blocks of equal angular size."};
};

struct CubeScale {
using type = double;
static constexpr Options::String help = {
"Specify the desired cube scale that must be greater than or equal to "
"1.0. The initial separation is multiplied by this cube scale to "
"produce larger cubes around each object which is desirable when "
"closer to merger."};
static double lower_bound() { return 1.0; }
};

struct InitialRefinement {
using type =
std::variant<size_t, std::array<size_t, 3>,
Expand Down Expand Up @@ -441,7 +459,7 @@ class BinaryCompactObject : public DomainCreator<3> {
template <typename Metavariables>
using options = tmpl::append<
tmpl::list<ObjectA, ObjectB, CenterOfMassOffset, EnvelopeRadius,
OuterRadius, InitialRefinement, InitialGridPoints,
OuterRadius, CubeScale, InitialRefinement, InitialGridPoints,
UseEquiangularMap, RadialDistributionEnvelope,
RadialDistributionOuterShell, OpeningAngle, TimeDependentMaps>,
tmpl::conditional_t<
Expand Down Expand Up @@ -480,7 +498,7 @@ class BinaryCompactObject : public DomainCreator<3> {
BinaryCompactObject(
typename ObjectA::type object_A, typename ObjectB::type object_B,
std::array<double, 2> center_of_mass_offset, double envelope_radius,
double outer_radius,
double outer_radius, double cube_scale,
const typename InitialRefinement::type& initial_refinement,
const typename InitialGridPoints::type& initial_number_of_grid_points,
bool use_equiangular_map = true,
Expand Down Expand Up @@ -560,6 +578,8 @@ class BinaryCompactObject : public DomainCreator<3> {
block_groups_{};
std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
grid_anchors_{};
double offset_x_coord_a_{};
double offset_x_coord_b_{};

// Variables to handle std::variant on Object A and B
double x_coord_a_{};
Expand Down
2 changes: 1 addition & 1 deletion src/Domain/Creators/CylindricalBinaryCompactObject.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -416,7 +416,7 @@ CylindricalBinaryCompactObject::CylindricalBinaryCompactObject(
time_dependent_options_->build_maps(
std::array{rotate_from_z_to_x_axis(center_A_),
rotate_from_z_to_x_axis(center_B_)},
std::array{radius_A_, outer_radius_A_},
std::nullopt, std::nullopt, std::array{radius_A_, outer_radius_A_},
std::array{radius_B_, outer_radius_B_}, inner_common_radius,
outer_radius_);
}
Expand Down
2 changes: 1 addition & 1 deletion src/Domain/Creators/Sphere.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -308,7 +308,7 @@ Domain<3> Sphere::create_domain() const {
sph_wedge_coordinate_maps(
inner_radius_, outer_radius_,
fill_interior_ ? std::get<InnerCube>(interior_).sphericity : 1.0, 1.0,
use_equiangular_map_, false, radial_partitioning_,
use_equiangular_map_, std::nullopt, false, radial_partitioning_,
radial_distribution_, which_wedges_),
compression);

Expand Down
Loading

0 comments on commit 66fb37d

Please sign in to comment.