Merge pull request #6212 from iago-mendes/control-adm

Control center of mass and linear momentum in initial data.

src/Domain/Creators/BinaryCompactObject.cpp

@@ -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"
@@ -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,
@@ -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),
@@ -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{
@@ -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_);
   }
 }
@@ -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) ---
     //
@@ -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) ---
     //
@@ -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_)
@@ -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_) {
@@ -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()}}});
     }
   }
 

src/Domain/Creators/BinaryCompactObject.hpp

@@ -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 =
@@ -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
@@ -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."};
@@ -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>,
@@ -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,
@@ -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_{};

src/PointwiseFunctions/AnalyticData/Xcts/Binary.hpp

@@ -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";
@@ -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.";
@@ -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),
@@ -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_;
@@ -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()
@@ -440,6 +455,8 @@ class Binary : public elliptic::analytic_data::Background,
 
  private:
   std::array<double, 2> xcoords_{};
+  double y_offset_{};
+  double z_offset_{};
   std::array<std::unique_ptr<IsolatedObjectBase>, 2> superposed_objects_{};
   Xcts::Solutions::Flatness flatness_{};
   double angular_velocity_ = std::numeric_limits<double>::signaling_NaN();
@@ -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
@@ -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)) /