WIP. new transition inside excision. Still needs dox

src/Domain/CoordinateMaps/TimeDependent/Shape.cpp

@@ -50,7 +50,7 @@ void Shape::jacobian_helper(
     gsl::not_null<tnsr::Ij<T, 3, Frame::NoFrame>*> result,
     const ylm::Spherepack::InterpolationInfo<T>& interpolation_info,
     const DataVector& extended_coefs, const std::array<T, 3>& centered_coords,
-    const T& radial_distortion, const T& transition_func_over_radius) const {
+    const T& radial_distortion, const T& transition_func) const {
   const auto angular_gradient =
       extended_ylm_.gradient_from_coefs(extended_coefs);
 
@@ -97,13 +97,13 @@ void Shape::jacobian_helper(
                               interpolation_info);
   }
 
-  auto transition_func_over_square_radius =
+  auto transition_func_over_radius =
       transition_func_->operator()(centered_coords, {1});
   auto transition_func_gradient_times_distortion =
       transition_func_->gradient(centered_coords) * radial_distortion;
 
   auto& target_gradient_times_spatial_part = target_gradient;
-  target_gradient_times_spatial_part *= transition_func_over_square_radius;
+  target_gradient_times_spatial_part *= transition_func_over_radius;
 
   for (size_t i = 0; i < 3; i++) {
     for (size_t j = 0; j < 3; j++) {
@@ -113,7 +113,7 @@ void Shape::jacobian_helper(
            gsl::at(target_gradient_times_spatial_part, j));
     }
 
-    result->get(i, i) += 1.0 - radial_distortion * transition_func_over_radius;
+    result->get(i, i) += 1.0 - radial_distortion * transition_func;
   }
 }
 
@@ -278,14 +278,13 @@ tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> Shape::jacobian(
                                        extended_coefs, interpolation_info);
 
   using ReturnType = tt::remove_cvref_wrap_t<T>;
-  const ReturnType transition_func_over_radius =
+  const ReturnType transition_func =
       transition_func_->operator()(centered_coords, std::nullopt);
   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> result(
       get_size(centered_coords[0]));
 
   jacobian_helper(make_not_null(&result), interpolation_info, extended_coefs,
-                  centered_coords, radial_distortion,
-                  transition_func_over_radius);
+                  centered_coords, radial_distortion, transition_func);
   return result;
 }
 
@@ -351,23 +350,20 @@ void Shape::coords_frame_velocity_jacobian(
   extended_ylm_.interpolate_from_coefs(make_not_null(&radial_distortion),
                                        extended_coefs, interpolation_info);
 
-  auto& transition_func_over_radius = get(get<::Tags::TempScalar<1>>(temps));
-  transition_func_over_radius =
-      transition_func_->operator()(centered_coords, std::nullopt);
+  auto& transition_func = get(get<::Tags::TempScalar<1>>(temps));
+  transition_func = transition_func_->operator()(centered_coords, std::nullopt);
   *source_and_target_coords =
-      center_ +
-      centered_coords * (1. - radial_distortion * transition_func_over_radius);
+      center_ + centered_coords * (1. - radial_distortion * transition_func);
 
   auto& radii_velocities = get<0, 1>(*jac);
   extended_ylm_.interpolate_from_coefs(make_not_null(&radii_velocities),
                                        extended_coefs_derivs,
                                        interpolation_info);
-  *frame_vel =
-      -centered_coords * radii_velocities * transition_func_over_radius;
+  *frame_vel = -centered_coords * radii_velocities * transition_func;
 
   jacobian_helper<DataVector>(jac, interpolation_info, extended_coefs,
                               centered_coords, radial_distortion,
-                              transition_func_over_radius);
+                              transition_func);
 }
 
 template <typename T>

src/Domain/CoordinateMaps/TimeDependent/Shape.hpp

@@ -300,7 +300,7 @@ class Shape {
       gsl::not_null<tnsr::Ij<T, 3, Frame::NoFrame>*> result,
       const ylm::Spherepack::InterpolationInfo<T>& interpolation_info,
       const DataVector& extended_coefs, const std::array<T, 3>& centered_coords,
-      const T& radial_distortion, const T& transition_func_over_radius) const;
+      const T& radial_distortion, const T& transition_func) const;
 
   void check_size(const gsl::not_null<DataVector*>& coefs,
                   const FunctionsOfTimeMap& functions_of_time, double time,

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.cpp

@@ -10,6 +10,7 @@
 
 #include "DataStructures/Blaze/IntegerPow.hpp"
 #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp"
+#include "Utilities/Algorithm.hpp"
 #include "Utilities/ConstantExpressions.hpp"
 #include "Utilities/ContainerHelpers.hpp"
 #include "Utilities/EqualWithinRoundoff.hpp"
@@ -27,6 +28,7 @@ SphereTransition::SphereTransition(const double r_min, const double r_max,
   if (r_min <= 0.) {
     ERROR("The minimum radius must be greater than 0 but is " << r_min);
   }
+  inverse_cube_r_min_ = 1.0 / cube(r_min_);
   if (r_max <= r_min) {
     ERROR(
         "The maximum radius must be greater than the minimum radius but "
@@ -49,17 +51,22 @@ double SphereTransition::operator()(
     const std::optional<size_t>& one_over_radius_power) const {
   const double mag = magnitude(source_coords);
 
-  if (UNLIKELY(one_over_radius_power.has_value() and
+  if (UNLIKELY(one_over_radius_power.value_or(0_st) >= 3 and
                equal_within_roundoff(mag, 0.0))) {
     ERROR("Trying to divide by a point "
           << source_coords
           << " with radius zero in SphereTransition operator.");
   }
 
   if (interior_) {
-    return 1.0 / (r_min_ *
-                  integer_pow(mag, static_cast<int>(
-                                       one_over_radius_power.value_or(0_st))));
+    return inverse_cube_r_min_ *
+           (one_over_radius_power.value_or(0_st) < 3
+                ? integer_pow(mag,
+                              static_cast<int>(
+                                  2 - one_over_radius_power.value_or(0_st)))
+                : 1.0 /
+                      integer_pow(mag, static_cast<int>(
+                                           one_over_radius_power.value() - 2)));
   }
 
   if (UNLIKELY(mag < r_min_ - eps_)) {
@@ -107,8 +114,45 @@ std::optional<double> SphereTransition::original_radius_over_radius(
   // If we aren't reversed, check near or within r_min_.
   if (a_ < 0.0) {
     if (mag + radial_distortion < r_min_ - eps_) {
-      return interior_ ? std::optional{r_min_ / (r_min_ - radial_distortion)}
-                       : std::nullopt;
+      if (not interior_) {
+        return std::nullopt;
+      }
+
+      const double factor =
+          1.0 / (inverse_cube_r_min_ * square(mag) * radial_distortion);
+
+      // The following variable names are defined in Numerical Recipes (3rd
+      // edition) pg 228. We choose to keep them as they appear in NR for better
+      // comparison with the book.
+
+      // Numerical Recipes (3rd edition) pg 228, eq 5.6.10
+      const double R = 0.5 * factor;
+      const double Q = factor / 3.0;
+      const bool multiple_real_roots = square(R) < cube(Q);
+
+      if (multiple_real_roots) {
+        // Numerical Recipes (3rd edition) pg 228, eqs 5.6.11 - 5.6.12
+        const double theta = acos(R / sqrt(cube(Q)));
+        std::vector<double> roots{
+            -2.0 * sqrt(Q) * cos(theta / 3.0),
+            -2.0 * sqrt(Q) * cos((theta + 2.0 * M_PI) / 3.0),
+            -2.0 * sqrt(Q) * cos((theta - 2.0 * M_PI) / 3.0)};
+
+        // Radii are positive
+        alg::remove_if(roots, [](const double root) { return root < 0.0; });
+
+        // Since the root of this is the original radius over target radius, it
+        // will be of order unity and not something super large, so we take the
+        // smallest of the positive roots. If this turns out to not be robust
+        // enough we can change this.
+        return *alg::min_element(roots);
+      } else {
+        // Numerical Recipes (3rd edition) pg 228, eqs 5.6.13 - 5.6.17
+        const double A = -sgn(R) * cbrt(abs(R) + sqrt(square(R) - cube(Q)));
+        const double B = A == 0.0 ? 0.0 : Q / A;
+
+        return A + B;
+      }
     } else if (equal_within_roundoff(mag, r_min_)) {
       return std::optional{r_min_ / (r_min_ - radial_distortion)};
     }
@@ -158,7 +202,7 @@ std::array<double, 3> SphereTransition::gradient(
 
   // Short circuit for the interior
   if (interior_) {
-    return std::array{0.0, 0.0, 0.0};
+    return 2.0 * inverse_cube_r_min_ * source_coords;
   }
 
   if (UNLIKELY(equal_within_roundoff(mag, 0.0))) {
@@ -225,11 +269,14 @@ void SphereTransition::pup(PUP::er& p) {
     p | r_max_;
     p | a_;
     p | b_;
-    if (version >= 1) {
-      p | interior_;
-    } else {
-      interior_ = false;
-    }
+  }
+
+  if (p.isUnpacking()) {
+    inverse_cube_r_min_ = 1.0 / cube(r_min_);
+  }
+
+  if (version >= 1) {
+    p | interior_;
   } else if (p.isUnpacking()) {
     interior_ = false;
   }

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.hpp

@@ -81,6 +81,7 @@ class SphereTransition final : public ShapeMapTransitionFunction {
 
  private:
   double r_min_{};
+  double inverse_cube_r_min_{};
   double r_max_{};
   double a_{};
   double b_{};

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/Wedge.cpp

@@ -339,7 +339,7 @@ double Wedge::operator()(
 
   using ::operator<<;
 
-  if (UNLIKELY(one_over_radius_power.has_value() and
+  if (UNLIKELY(one_over_radius_power.value_or(0_st) >= 3 and
                equal_within_roundoff(centered_coords_magnitude, 0.0))) {
     ERROR("Trying to divide by a point "
           << source_coords
@@ -348,10 +348,14 @@ double Wedge::operator()(
 
   // If the axis is the interior, short circuit
   if (axis_ == Axis::Interior) {
-    return 1.0 / (inner_surface_.radius *
-                  integer_pow(
-                      centered_coords_magnitude,
-                      static_cast<int>(one_over_radius_power.value_or(0_st))));
+    return 1.0 / cube(inner_surface_.radius) *
+           (one_over_radius_power.value_or(0_st) < 3
+                ? integer_pow(centered_coords_magnitude,
+                              static_cast<int>(
+                                  2 - one_over_radius_power.value_or(0_st)))
+                : 1.0 / integer_pow(centered_coords_magnitude,
+                                    static_cast<int>(
+                                        one_over_radius_power.value() - 2)));
   }
 
   // First check if we are at the inner center to avoid dividing by zero below
@@ -451,14 +455,50 @@ std::optional<double> Wedge::original_radius_over_radius(
 
   if (not reverse_) {
     // Check if we are within the interior region or at the inner surface. The
-    // formula is simplified in this case.
+    // formula is simplified at the inner surface
     if (centered_coords_magnitude + radial_distortion < inner_distance - eps_) {
       // Int: Unless the axis_ is Interior, this block doesn't contain this
       // point
-      return axis_ == Axis::Interior
-                 ? std::optional{inner_distance /
-                                 (inner_distance - radial_distortion)}
-                 : std::nullopt;
+      if (axis_ != Axis::Interior) {
+        return std::nullopt;
+      }
+
+      const double factor =
+          cube(inner_surface_.radius) /
+          (square(centered_coords_magnitude) * radial_distortion);
+
+      // The following variable names are defined in Numerical Recipes (3rd
+      // edition) pg 228. We choose to keep them as they appear in NR for better
+      // comparison with the book.
+
+      // Numerical Recipes (3rd edition) pg 228, eq 5.6.10
+      const double R = 0.5 * factor;
+      const double Q = factor / 3.0;
+      const bool multiple_real_roots = square(R) < cube(Q);
+
+      if (multiple_real_roots) {
+        // Numerical Recipes (3rd edition) pg 228, eqs 5.6.11 - 5.6.12
+        const double theta = acos(R / sqrt(cube(Q)));
+        std::vector<double> roots{
+            -2.0 * sqrt(Q) * cos(theta / 3.0),
+            -2.0 * sqrt(Q) * cos((theta + 2.0 * M_PI) / 3.0),
+            -2.0 * sqrt(Q) * cos((theta - 2.0 * M_PI) / 3.0)};
+
+        // Radii are positive
+        alg::remove_if(roots, [](const double root) { return root < 0.0; });
+
+        // Since the root of this is the original radius over target radius, it
+        // will be of order unity and not something super large, so we take the
+        // smallest of the positive roots. If this turns out to not be robust
+        // enough we can change this.
+        return *alg::min_element(roots);
+      } else {
+        // Numerical Recipes (3rd edition) pg 228, eqs 5.6.13 - 5.6.17
+        const double A = -sgn(R) * cbrt(abs(R) + sqrt(square(R) - cube(Q)));
+        const double B = A == 0.0 ? 0.0 : Q / A;
+
+        return A + B;
+      }
     } else if (equal_within_roundoff(
                    centered_coords_magnitude + radial_distortion,
                    inner_distance)) {
@@ -539,10 +579,9 @@ std::array<double, 3> Wedge::gradient(
   // The source coords are centered
   const double centered_coords_magnitude = magnitude(source_coords);
 
-  // Short circuit if we are in the interior and just return zero for the
-  // gradient
+  // Short circuit if we are in the interior
   if (axis_ == Axis::Interior) {
-    return std::array{0.0, 0.0, 0.0};
+    return 2.0 / cube(inner_surface_.radius) * source_coords;
   }
 
   // First check if we are at the inner center to avoid dividing by zero below

tests/Unit/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/Test_SphereTransition.cpp

@@ -1,6 +1,7 @@
 // Distributed under the MIT License.
 // See LICENSE.txt for details.
 
+#include "Framework/TestHelpers.hpp"
 #include "Framework/TestingFramework.hpp"
 
 #include <array>
@@ -39,8 +40,11 @@ SPECTRE_TEST_CASE("Unit.Domain.CoordinateMaps.Shape.SphereTransition",
 
   {
     INFO("Sphere transition");
-    const SphereTransition sphere_transition{2., 4.};
-    const SphereTransition sphere_transition_interior{2., 4., false, true};
+    SphereTransition sphere_transition{2., 4.};
+    sphere_transition = serialize_and_deserialize(sphere_transition);
+    SphereTransition sphere_transition_interior{2., 4., false, true};
+    sphere_transition_interior =
+        serialize_and_deserialize(sphere_transition_interior);
     const domain::CoordinateMaps::TimeDependent::Shape shape_map{
         std::array{0.0, 0.0, 0.0}, l_max, l_max,
         std::make_unique<SphereTransition>(sphere_transition), "Shape"};
@@ -64,16 +68,25 @@ SPECTRE_TEST_CASE("Unit.Domain.CoordinateMaps.Shape.SphereTransition",
 
     const std::vector<std::array<double, 3>> interior_points{{1.0, 0.0, 0.0},
                                                              {0.0, 0.0, 0.0}};
-    for (const auto& interior_point : interior_points) {
-      CAPTURE(interior_point);
-      CHECK(sphere_transition_interior(interior_point, std::nullopt) == 0.5);
-      test_inverse_map(shape_map_interior, interior_point, time,
+    const std::vector<double> interior_function_values{0.125, 0.0};
+    for (size_t i = 0; i < interior_points.size(); i++) {
+      CAPTURE(interior_points[i]);
+      CHECK(sphere_transition_interior(interior_points[i], std::nullopt) ==
+            interior_function_values[i]);
+      test_inverse_map(shape_map_interior, interior_points[i], time,
                        functions_of_time);
     }
+
+    // Check close, but not at, center and r_min
+    test_inverse_map(shape_map_interior, std::array{2.0 * eps, 0.0, 0.0}, time,
+                     functions_of_time);
+    test_inverse_map(shape_map_interior, std::array{2.0 - 2.0 * eps, 0.0, 0.0},
+                     time, functions_of_time);
   }
   {
     INFO("Reverse sphere transition");
-    const SphereTransition sphere_transition{2., 4., true};
+    SphereTransition sphere_transition{2., 4., true};
+    sphere_transition = serialize_and_deserialize(sphere_transition);
 
     const domain::CoordinateMaps::TimeDependent::Shape shape_map{
         std::array{0.0, 0.0, 0.0}, 4, 4,