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

Control center of mass and linear momentum in initial data. #6212

Merged
merged 2 commits into from
Aug 16, 2024
Merged

Control center of mass and linear momentum in initial data. #6212

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b865cbd
Prepare BBH domain for center-of-mass control.
iago-mendes Aug 15, 2024
ee7e100
Control center of mass and linear momentum in initial data.
iago-mendes Aug 15, 2024
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
107 changes: 64 additions & 43 deletions src/Domain/Creators/BinaryCompactObject.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,6 @@
#include "Domain/CoordinateMaps/Distribution.hpp"
#include "Domain/CoordinateMaps/Equiangular.hpp"
#include "Domain/CoordinateMaps/Frustum.hpp"
#include "Domain/CoordinateMaps/Identity.hpp"
#include "Domain/CoordinateMaps/Interval.hpp"
#include "Domain/CoordinateMaps/ProductMaps.hpp"
#include "Domain/CoordinateMaps/ProductMaps.tpp"
Expand Down Expand Up @@ -71,7 +70,8 @@ bool BinaryCompactObject::Object::is_excised() const {

BinaryCompactObject::BinaryCompactObject(
typename ObjectA::type object_A, typename ObjectB::type object_B,
const double envelope_radius, const double outer_radius,
std::array<double, 2> center_of_mass_offset, const double envelope_radius,
const double outer_radius,
const typename InitialRefinement::type& initial_refinement,
const typename InitialGridPoints::type& initial_number_of_grid_points,
const bool use_equiangular_map,
Expand All @@ -84,6 +84,7 @@ BinaryCompactObject::BinaryCompactObject(
const Options::Context& context)
: object_A_(std::move(object_A)),
object_B_(std::move(object_B)),
center_of_mass_offset_(center_of_mass_offset),
envelope_radius_(envelope_radius),
outer_radius_(outer_radius),
use_equiangular_map_(use_equiangular_map),
Expand Down Expand Up @@ -215,8 +216,11 @@ BinaryCompactObject::BinaryCompactObject(
}

// Create grid anchors
grid_anchors_ = bco::create_grid_anchors(std::array{x_coord_a_, 0.0, 0.0},
std::array{x_coord_b_, 0.0, 0.0});
grid_anchors_ =
bco::create_grid_anchors(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]});

// Create block names and groups
static std::array<std::string, 6> wedge_directions{
Expand Down Expand Up @@ -324,8 +328,10 @@ BinaryCompactObject::BinaryCompactObject(

if (time_dependent_options_.has_value()) {
time_dependent_options_->build_maps(
std::array{std::array{x_coord_a_, 0.0, 0.0},
std::array{x_coord_b_, 0.0, 0.0}},
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_);
}
}
Expand Down Expand Up @@ -354,21 +360,31 @@ Domain<3> BinaryCompactObject::create_domain() const {
Maps maps{};

// ObjectA/B is on the right/left, respectively.
const Translation translation_A{
Affine{-1.0, 1.0, -1.0 + x_coord_a_, 1.0 + x_coord_a_}, Identity2D{}};
const Translation translation_B{
Affine{-1.0, 1.0, -1.0 + x_coord_b_, 1.0 + x_coord_b_}, Identity2D{}};
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 + 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 + 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]}};

// Two blocks covering the compact objects and their immediate neighborhood
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_,
0.5 * length_inner_cube_),
Affine(-1.0, 1.0, -0.5 * length_inner_cube_,
0.5 * length_inner_cube_)}));
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_ + center_of_mass_offset_[0],
0.5 * length_inner_cube_ + center_of_mass_offset_[0]),
Affine(-1.0, 1.0,
-0.5 * length_inner_cube_ + center_of_mass_offset_[1],
0.5 * length_inner_cube_ + center_of_mass_offset_[1])}));
} else {
// --- Blocks enclosing each object (12 blocks per object) ---
//
Expand Down Expand Up @@ -397,13 +413,15 @@ Domain<3> BinaryCompactObject::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_,
0.5 * length_inner_cube_),
Affine(-1.0, 1.0, -0.5 * length_inner_cube_,
0.5 * length_inner_cube_)}));
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_ + center_of_mass_offset_[0],
0.5 * length_inner_cube_ + center_of_mass_offset_[0]),
Affine(-1.0, 1.0,
-0.5 * length_inner_cube_ + center_of_mass_offset_[1],
0.5 * length_inner_cube_ + center_of_mass_offset_[1])}));
} else {
// --- Blocks enclosing each object (12 blocks per object) ---
//
Expand Down Expand Up @@ -439,7 +457,8 @@ Domain<3> BinaryCompactObject::create_domain() const {
Maps maps_frustums = domain::make_vector_coordinate_map_base<
Frame::BlockLogical, Frame::Inertial, 3>(frustum_coordinate_maps(
length_inner_cube_, length_outer_cube_, use_equiangular_map_,
{{-translation_, 0.0, 0.0}}, radial_distribution_envelope_,
{{-translation_, -center_of_mass_offset_[0], -center_of_mass_offset_[1]}},
radial_distribution_envelope_,
radial_distribution_envelope_ ==
domain::CoordinateMaps::Distribution::Projective
? std::optional<double>(length_inner_cube_ / length_outer_cube_)
Expand Down Expand Up @@ -494,15 +513,16 @@ Domain<3> BinaryCompactObject::create_domain() const {
else {
excision_spheres.emplace(
"ExcisionSphereA",
ExcisionSphere<3>{
std::get<Object>(object_A_).inner_radius,
tnsr::I<double, 3, Frame::Grid>{{x_coord_a_, 0.0, 0.0}},
{{0, Direction<3>::lower_zeta()},
{1, Direction<3>::lower_zeta()},
{2, Direction<3>::lower_zeta()},
{3, Direction<3>::lower_zeta()},
{4, Direction<3>::lower_zeta()},
{5, Direction<3>::lower_zeta()}}});
ExcisionSphere<3>{std::get<Object>(object_A_).inner_radius,
tnsr::I<double, 3, Frame::Grid>{
{x_coord_a_, center_of_mass_offset_[0],
center_of_mass_offset_[1]}},
{{0, Direction<3>::lower_zeta()},
{1, Direction<3>::lower_zeta()},
{2, Direction<3>::lower_zeta()},
{3, Direction<3>::lower_zeta()},
{4, Direction<3>::lower_zeta()},
{5, Direction<3>::lower_zeta()}}});
}
}
if (not use_single_block_b_) {
Expand Down Expand Up @@ -536,15 +556,16 @@ Domain<3> BinaryCompactObject::create_domain() const {
else {
excision_spheres.emplace(
"ExcisionSphereB",
ExcisionSphere<3>{
std::get<Object>(object_B_).inner_radius,
tnsr::I<double, 3, Frame::Grid>{{x_coord_b_, 0.0, 0.0}},
{{12, Direction<3>::lower_zeta()},
{13, Direction<3>::lower_zeta()},
{14, Direction<3>::lower_zeta()},
{15, Direction<3>::lower_zeta()},
{16, Direction<3>::lower_zeta()},
{17, Direction<3>::lower_zeta()}}});
ExcisionSphere<3>{std::get<Object>(object_B_).inner_radius,
tnsr::I<double, 3, Frame::Grid>{
{x_coord_b_, center_of_mass_offset_[0],
center_of_mass_offset_[1]}},
{{12, Direction<3>::lower_zeta()},
{13, Direction<3>::lower_zeta()},
{14, Direction<3>::lower_zeta()},
{15, Direction<3>::lower_zeta()},
{16, Direction<3>::lower_zeta()},
{17, Direction<3>::lower_zeta()}}});
}
}

Expand Down
29 changes: 24 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 @@ -151,6 +151,8 @@ class BinaryCompactObject : public DomainCreator<3> {
using Affine = CoordinateMaps::Affine;
using Affine3D = CoordinateMaps::ProductOf3Maps<Affine, Affine, Affine>;
using Identity2D = CoordinateMaps::Identity<2>;
// The Translation type is no longer needed, but it is kept here for backwards
// compatibility with old domains.
using Translation = CoordinateMaps::ProductOf2Maps<Affine, Identity2D>;
using Equiangular = CoordinateMaps::Equiangular;
using Equiangular3D =
Expand All @@ -175,6 +177,12 @@ class BinaryCompactObject : public DomainCreator<3> {
CoordinateMaps::Wedge<3>>,
domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
CoordinateMaps::Wedge<3>, Translation>,
domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Affine3D,
Affine3D>,
domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Equiangular3D,
Affine3D>,
domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
CoordinateMaps::Wedge<3>, Affine3D>,
bco::TimeDependentMapOptions<false>::maps_list>>;

/// Options for an excision region in the domain
Expand Down Expand Up @@ -311,6 +319,15 @@ class BinaryCompactObject : public DomainCreator<3> {
"x-axis)."};
};

struct CenterOfMassOffset {
using type = std::array<double, 2>;
static constexpr Options::String help = {
"Offset in the y and z axes applied to both object A and B in order to "
"control the center of mass. This moves the location of the two objects"
" in the grid frame but keeps the Envelope and OuterShell centered on "
"the origin in the grid frame."};
};

struct Envelope {
static constexpr Options::String help = {
"Options for the sphere enveloping the two objects."};
Expand Down Expand Up @@ -410,10 +427,10 @@ class BinaryCompactObject : public DomainCreator<3> {

template <typename Metavariables>
using options = tmpl::append<
tmpl::list<ObjectA, ObjectB, EnvelopeRadius, OuterRadius,
InitialRefinement, InitialGridPoints, UseEquiangularMap,
RadialDistributionEnvelope, RadialDistributionOuterShell,
OpeningAngle, TimeDependentMaps>,
tmpl::list<ObjectA, ObjectB, CenterOfMassOffset, EnvelopeRadius,
OuterRadius, InitialRefinement, InitialGridPoints,
UseEquiangularMap, RadialDistributionEnvelope,
RadialDistributionOuterShell, OpeningAngle, TimeDependentMaps>,
tmpl::conditional_t<
domain::BoundaryConditions::has_boundary_conditions_base_v<
typename Metavariables::system>,
Expand Down Expand Up @@ -449,7 +466,8 @@ class BinaryCompactObject : public DomainCreator<3> {

BinaryCompactObject(
typename ObjectA::type object_A, typename ObjectB::type object_B,
double envelope_radius, double outer_radius,
std::array<double, 2> center_of_mass_offset, double envelope_radius,
double outer_radius,
const typename InitialRefinement::type& initial_refinement,
const typename InitialGridPoints::type& initial_number_of_grid_points,
bool use_equiangular_map = true,
Expand Down Expand Up @@ -507,6 +525,7 @@ class BinaryCompactObject : public DomainCreator<3> {
private:
typename ObjectA::type object_A_{};
typename ObjectB::type object_B_{};
std::array<double, 2> center_of_mass_offset_{};
double envelope_radius_ = std::numeric_limits<double>::signaling_NaN();
double outer_radius_ = std::numeric_limits<double>::signaling_NaN();
std::vector<std::array<size_t, 3>> initial_refinement_{};
Expand Down
28 changes: 25 additions & 3 deletions src/PointwiseFunctions/AnalyticData/Xcts/Binary.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -321,6 +321,12 @@ class Binary : public elliptic::analytic_data::Background,
"The coordinates on the x-axis where the two objects are placed";
using type = std::array<double, 2>;
};
struct CenterOfMassOffset {
static constexpr Options::String help = {
"Offset in the y and z axes applied to both objects in order to "
"control the center of mass."};
using type = std::array<double, 2>;
};
struct ObjectLeft {
static constexpr Options::String help =
"The object placed on the negative x-axis";
Expand Down Expand Up @@ -355,8 +361,9 @@ class Binary : public elliptic::analytic_data::Background,
"to disable the Gaussian falloff.";
using type = Options::Auto<std::array<double, 2>, Options::AutoLabel::None>;
};
using options = tmpl::list<XCoords, ObjectLeft, ObjectRight, AngularVelocity,
Expansion, LinearVelocity, FalloffWidths>;
using options =
tmpl::list<XCoords, CenterOfMassOffset, ObjectLeft, ObjectRight,
AngularVelocity, Expansion, LinearVelocity, FalloffWidths>;
static constexpr Options::String help =
"Binary compact-object data in general relativity, constructed from "
"superpositions of two isolated objects.";
Expand All @@ -369,13 +376,16 @@ class Binary : public elliptic::analytic_data::Background,
~Binary() = default;

Binary(const std::array<double, 2> xcoords,
const std::array<double, 2> center_of_mass_offset,
std::unique_ptr<IsolatedObjectBase> object_left,
std::unique_ptr<IsolatedObjectBase> object_right,
const double angular_velocity, const double expansion,
const std::array<double, 3> linear_velocity,
const std::optional<std::array<double, 2>> falloff_widths,
const Options::Context& context = {})
: xcoords_(xcoords),
y_offset_(center_of_mass_offset[0]),
z_offset_(center_of_mass_offset[1]),
superposed_objects_({std::move(object_left), std::move(object_right)}),
angular_velocity_(angular_velocity),
expansion_(expansion),
Expand Down Expand Up @@ -414,6 +424,8 @@ class Binary : public elliptic::analytic_data::Background,
elliptic::analytic_data::Background::pup(p);
elliptic::analytic_data::InitialGuess::pup(p);
p | xcoords_;
p | y_offset_;
p | z_offset_;
p | superposed_objects_;
p | angular_velocity_;
p | expansion_;
Expand All @@ -423,6 +435,9 @@ class Binary : public elliptic::analytic_data::Background,

/// Coordinates of the objects, ascending left to right
const std::array<double, 2>& x_coords() const { return xcoords_; }
/// Offset in y and z coordinates of the objects
double y_offset() const { return y_offset_; }
double z_offset() const { return z_offset_; }
/// The two objects. First entry is the left object, second entry is the right
/// object.
const std::array<std::unique_ptr<IsolatedObjectBase>, 2>& superposed_objects()
Expand All @@ -440,6 +455,8 @@ class Binary : public elliptic::analytic_data::Background,

private:
std::array<double, 2> xcoords_{};
double y_offset_{};
double z_offset_{};
Comment on lines +458 to +459
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These are probably ok as is

std::array<std::unique_ptr<IsolatedObjectBase>, 2> superposed_objects_{};
Xcts::Solutions::Flatness flatness_{};
double angular_velocity_ = std::numeric_limits<double>::signaling_NaN();
Expand All @@ -456,6 +473,9 @@ class Binary : public elliptic::analytic_data::Background,
inv_jacobian,
tmpl::list<RequestedTags...> /*meta*/) const {
std::array<tnsr::I<DataVector, 3>, 2> x_isolated{{x, x}};
const std::array<std::array<double, 3>, 2> coords_isolated{
{{{xcoords_[0], y_offset_, z_offset_}},
{{xcoords_[1], y_offset_, z_offset_}}}};
std::array<DataVector, 2> euclidean_distance{};
std::array<DataVector, 2> windows{};
// Possible optimization: Only retrieve those superposed tags from the
Expand All @@ -466,7 +486,9 @@ class Binary : public elliptic::analytic_data::Background,
std::array<tuples::tagged_tuple_from_typelist<requested_superposed_tags>, 2>
isolated_vars;
for (size_t i = 0; i < 2; ++i) {
get<0>(gsl::at(x_isolated, i)) -= gsl::at(xcoords_, i);
for (size_t dim = 0; dim < 3; dim++) {
gsl::at(x_isolated, i).get(dim) -= gsl::at(coords_isolated, i)[dim];
}
gsl::at(euclidean_distance, i) = get(magnitude(gsl::at(x_isolated, i)));
if (falloff_widths_.has_value()) {
gsl::at(windows, i) = exp(-square(gsl::at(euclidean_distance, i)) /
Expand Down
Loading
Loading