Domains with large outer radius are inconsistent

[nilsvu]

Bug reports:

Expected behavior:

I can create a Sphere or BinaryCompactObject domain with an outer radius of 1e11.

Current behavior:

Relevant code: Wedge coordinate map https://github.com/sxs-collaboration/spectre/blob/develop/src/Domain/CoordinateMaps/Wedge.hpp

The following test fails:

diff --git a/tests/Unit/Domain/Creators/Test_Sphere.cpp b/tests/Unit/Domain/Creators/Test_Sphere.cpp
index 2db322e174..f538623f06 100644
--- a/tests/Unit/Domain/Creators/Test_Sphere.cpp
+++ b/tests/Unit/Domain/Creators/Test_Sphere.cpp
@@ -679,5 +679,19 @@ SPECTRE_TEST_CASE("Unit.Domain.Creators.Sphere", "[Domain][Unit]") {
   domain::creators::time_dependence::register_derived_with_charm();
   test_parse_errors();
   test_sphere();
+  const creators::Sphere sphere{
+      1.,
+      1e11,  // large outer radius!
+      creators::Sphere::Excision{create_boundary_condition(true)},
+      0_st,
+      8_st,
+      true,
+      std::nullopt,
+      {},
+      CoordinateMaps::Distribution::Inverse,
+      ShellWedges::All,
+      std::nullopt,
+      create_boundary_condition(true)};
+  TestHelpers::domain::creators::test_domain_creator(sphere, true);
 }
 }  // namespace domain

The test complains that the domain::physical_separation is too large. That might just be an issue with tolerances. However, the bigger issue is that tolerances like this are also used in production code when determining corners. Make the following change:

diff --git a/src/Domain/Creators/Sphere.cpp b/src/Domain/Creators/Sphere.cpp
index b244b51f24..d68790fea5 100644
--- a/src/Domain/Creators/Sphere.cpp
+++ b/src/Domain/Creators/Sphere.cpp
@@ -350,7 +350,7 @@ Domain<3> Sphere::create_domain() const {
     }
   }

-  Domain<3> domain{std::move(coord_maps),       corners,      {},
+  Domain<3> domain{std::move(coord_maps),
                    std::move(excision_spheres), block_names_, block_groups_};
   ASSERT(domain.blocks().size() == num_blocks_,
          "Unexpected number of blocks. Expected "

Now the maps are used to determine block neighbors (like in the BBH domain). The test above will now fail with an incorrect number of boundary conditions. That's pretty catastrophic for the simulation :)

Finally, try this change:

diff --git a/tests/Unit/Domain/CoordinateMaps/Test_Wedge3D.cpp b/tests/Unit/Domain/CoordinateMaps/Test_Wedge3D.cpp
index 08dde17c48..c9428c5ba1 100644
--- a/tests/Unit/Domain/CoordinateMaps/Test_Wedge3D.cpp
+++ b/tests/Unit/Domain/CoordinateMaps/Test_Wedge3D.cpp
@@ -33,7 +33,7 @@ void test_wedge3d_all_directions() {
   std::uniform_real_distribution<> angle_dis(80.0, 90.0);
   const double inner_radius = inner_dis(gen);
   CAPTURE(inner_radius);
-  const double outer_radius = outer_dis(gen);
+  const double outer_radius = 1e11;
   CAPTURE(outer_radius);
   const double inner_sphericity = unit_dis(gen);
   CAPTURE(inner_sphericity);
(

Now the Wedge3D test fails with very wrong Jacobians. In particular, when the Inverse radial distribution is selected (note: the choice of radial distribution is randomized in the test). This looks like an issue with numerical cancellations, and it improves a bit with the following change:

diff --git a/src/Domain/CoordinateMaps/Wedge.cpp b/src/Domain/CoordinateMaps/Wedge.cpp
index 9ab76f5d7b..46dc91e398 100644
--- a/src/Domain/CoordinateMaps/Wedge.cpp
+++ b/src/Domain/CoordinateMaps/Wedge.cpp
@@ -113,7 +113,8 @@ tt::remove_cvref_wrap_t<T> Wedge<Dim>::default_physical_z(
   } else if (radial_distribution_ == Distribution::Logarithmic) {
     return exp(sphere_zero_ + sphere_rate_ * zeta) * one_over_rho;
   } else {
-    return one_over_rho / (sphere_zero_ + sphere_rate_ * zeta);
+    return 2. * one_over_rho /
+           ((1. + zeta) / radius_outer_ + (1. - zeta) / radius_inner_);
   }
 }

[nilsvu]

@mar-celine would be willing to help me fix this?

[nilsvu]

Partly addressed in #5606.

[knelli2]

@nilsvu Is this still a problem since you fixed ID with a large outer boundary?