Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add Wedge Offset Option to BCO #6243

Merged
merged 1 commit into from
Sep 7, 2024
Merged

Add Wedge Offset Option to BCO #6243

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file modified docs/Images/Domain/DomainCreators/binary_compact_object_domain.png
View file
Open in desktop
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
Open in desktop
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_)) {
knelli2 marked this conversation as resolved.
Show resolved Hide resolved
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading