Extend Sphere transition beyond boundaries

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

@@ -12,12 +12,13 @@
 #include "Utilities/ContainerHelpers.hpp"
 #include "Utilities/EqualWithinRoundoff.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
+#include "Utilities/MakeString.hpp"
 
 namespace domain::CoordinateMaps::ShapeMapTransitionFunctions {
 
 SphereTransition::SphereTransition(const double r_min, const double r_max,
-                                   const bool reverse)
-    : r_min_(r_min), r_max_(r_max) {
+                                   const bool reverse, const bool extend)
+    : r_min_(r_min), r_max_(r_max), extend_(extend) {
   if (r_min <= 0.) {
     ERROR("The minimum radius must be greater than 0 but is " << r_min);
   }
@@ -55,10 +56,26 @@ std::optional<double> SphereTransition::original_radius_over_radius(
   }
   const double original_radius = (mag + radial_distortion * b_) / denom;
 
-  return (original_radius + eps_) >= r_min_ and
-                 (original_radius - eps_) <= r_max_
-             ? std::optional<double>{original_radius / mag}
-             : std::nullopt;
+  // If we are within r_min and r_max, doesn't matter if we are reverse
+  // or extending, we just return the value we calculated. If we are extending,
+  // then if we are reverse and the point is outside r_max or we aren't reversed
+  // and the point is inside r_min, then we return a simplified formula since
+  // the transition function is 1 in this region. If we are extending, but the
+  // above conditions aren't true, then we are in the region where the
+  // transition function is 0, so we return 1.0. Otherwise we return nullopt.
+  if ((original_radius + eps_) >= r_min_ and
+      (original_radius - eps_) <= r_max_) {
+    return std::optional<double>{original_radius / mag};
+  } else if (extend_) {
+    // a_ being positive is a sentinel for reverse (see constructor)
+    if ((a_ > 0.0 and mag > r_max_) or (a_ < 0.0 and mag < r_min_)) {
+      return {1.0 + radial_distortion / mag};
+    } else {
+      return {1.0};
+    }
+  } else {
+    return std::nullopt;
+  }
 }
 
 std::array<double, 3> SphereTransition::gradient(
@@ -72,17 +89,29 @@ std::array<DataVector, 3> SphereTransition::gradient(
 
 template <typename T>
 T SphereTransition::call_impl(const std::array<T, 3>& source_coords) const {
-  const T mag = magnitude(source_coords);
+  T mag = magnitude(source_coords);
   check_magnitudes(mag);
-  return a_ * mag + b_;
+  auto& result = mag;
+  result = a_ * mag + b_;
+  return extend_ ? blaze::clamp(result, 0.0, 1.0) : result;
 }
 
 template <typename T>
 std::array<T, 3> SphereTransition::gradient_impl(
     const std::array<T, 3>& source_coords) const {
   const T mag = magnitude(source_coords);
   check_magnitudes(mag);
-  return a_ * source_coords / mag;
+  std::array<T, 3> result = a_ * source_coords / mag;
+  if (extend_) {
+    for (size_t i = 0; i < get_size(mag); i++) {
+      if (get_element(mag, i) < r_min_ or get_element(mag, i) > r_max_) {
+        for (size_t j = 0; j < 3; j++) {
+          get_element(gsl::at(result, j), i) = 0.0;
+        }
+      }
+    }
+  }
+  return result;
 }
 
 bool SphereTransition::operator==(
@@ -93,43 +122,65 @@ bool SphereTransition::operator==(
   const auto& derived = dynamic_cast<const SphereTransition&>(other);
   // no need to check `a_` and `b_` as they are uniquely determined by
   // `r_min_` and `r_max_`.
-  return this->r_min_ == derived.r_min_ and this->r_max_ == derived.r_max_;
+  return this->r_min_ == derived.r_min_ and this->r_max_ == derived.r_max_ and
+         this->extend_ == derived.extend_;
 }
 
 bool SphereTransition::operator!=(
     const ShapeMapTransitionFunction& other) const {
   return not(*this == other);
 }
 
-// checks that the magnitudes are all between `r_min_` and `r_max_`
+// checks that the magnitudes are all between `r_min_` and `r_max_` if we aren't
+// extending
 template <typename T>
 void SphereTransition::check_magnitudes([[maybe_unused]] const T& mag) const {
 #ifdef SPECTRE_DEBUG
   for (size_t i = 0; i < get_size(mag); ++i) {
-    if (get_element(mag, i) + eps_ < r_min_ or
-        get_element(mag, i) - eps_ > r_max_) {
+    const bool point_is_bad = extend_ ? get_element(mag, i) <= 0.0
+                                      : (get_element(mag, i) + eps_ < r_min_ or
+                                         get_element(mag, i) - eps_ > r_max_);
+    if (point_is_bad) {
       ERROR(
-          "The sphere transition map was called with coordinates outside the "
-          "set minimum and maximum radius. The minimum radius is "
-          << r_min_ << ", the maximum radius is " << r_max_
-          << ". The requested point has magnitude: " << get_element(mag, i));
+          "The sphere transition map was called with bad coordinates. The "
+          "requested point has magnitude "
+          << get_element(mag, i)
+          << (extend_
+                  ? (MakeString{} << " <= 0.0")
+                  : (MakeString{} << " which is outside the set minimum and "
+                                     "maxiumum radius. The minimum radius is "
+                                  << r_min_ << ", the maximum radius is "
+                                  << r_max_ << ".")));
     }
   }
 #endif  // SPECTRE_DEBUG
 }
 
 void SphereTransition::pup(PUP::er& p) {
   ShapeMapTransitionFunction::pup(p);
-  size_t version = 0;
+  size_t version = 1;
   p | version;
   // Remember to increment the version number when making changes to this
   // function. Retain support for unpacking data written by previous versions
   // whenever possible. See `Domain` docs for details.
-  if (version >= 0) {
+  if (p.isUnpacking()) {
+    if (version >= 0) {
+      p | r_min_;
+      p | r_max_;
+      p | a_;
+      p | b_;
+    }
+    if (version >= 1) {
+      p | extend_;
+    } else {
+      extend_ = false;
+    }
+  } else {
     p | r_min_;
     p | r_max_;
     p | a_;
     p | b_;
+    p | extend_;
   }
 }
 

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

@@ -26,14 +26,19 @@ namespace domain::CoordinateMaps::ShapeMapTransitionFunctions {
  * r_{\text{max}}
  * \f}
  *
- * If the `reverse` flag is set to `true`, then the function falls off from 0 at
+ * If \p reverse is set to `true`, then the function falls off from 0 at
  * `r_min` to 1 at `r_max`. To do this, the coefficients are modified as
  * $a \rightarrow -a$ and $b \rightarrow 1-b$.
+ *
+ * If \p extend is set to `true`, then the function can be called beyond `r_min`
+ * and `r_max`. Within `r_min` the value is 1, and outside `r_max` the value is
+ * 0. This is reversed if \p reverse is true.
  */
 class SphereTransition final : public ShapeMapTransitionFunction {
  public:
   explicit SphereTransition() = default;
-  SphereTransition(double r_min, double r_max, bool reverse = false);
+  SphereTransition(double r_min, double r_max, bool reverse = false,
+                   bool extend = false);
 
   double operator()(const std::array<double, 3>& source_coords) const override;
   DataVector operator()(
@@ -74,6 +79,7 @@ class SphereTransition final : public ShapeMapTransitionFunction {
   double r_max_{};
   double a_{};
   double b_{};
+  bool extend_{};
   static constexpr double eps_ = std::numeric_limits<double>::epsilon() * 100;
 };
 }  // namespace domain::CoordinateMaps::ShapeMapTransitionFunctions

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

@@ -26,6 +26,19 @@ SPECTRE_TEST_CASE("Unit.Domain.CoordinateMaps.Shape.SphereTransition",
     CHECK(sphere_transition(upper_bound) == approx(0.));
     const std::array<double, 3> upper_bound_eps{{4. + eps, 0., 0.}};
     CHECK(sphere_transition(upper_bound_eps) == approx(0.));
+    {
+      INFO("Extended");
+      const SphereTransition sphere_transition_extended{2., 4., false, true};
+      CHECK(sphere_transition_extended(lower_bound) == approx(1.0));
+      CHECK(sphere_transition_extended(lower_bound_eps) == approx(1.0));
+      const std::array<double, 3> inside_lower_bound{{1., 0., 0.}};
+      CHECK(sphere_transition_extended(inside_lower_bound) == 1.0);
+      CHECK(sphere_transition_extended(midpoint) == approx(0.5));
+      CHECK(sphere_transition_extended(upper_bound) == approx(0.));
+      CHECK(sphere_transition_extended(upper_bound_eps) == approx(0.));
+      const std::array<double, 3> outside_upper_bound{{5., 0., 0.}};
+      CHECK(sphere_transition_extended(outside_upper_bound) == 0.0);
+    }
   }
   {
     INFO("Reverse sphere transition");
@@ -40,6 +53,19 @@ SPECTRE_TEST_CASE("Unit.Domain.CoordinateMaps.Shape.SphereTransition",
     CHECK(sphere_transition(upper_bound) == approx(1.0));
     const std::array<double, 3> upper_bound_eps{{4. + eps, 0., 0.}};
     CHECK(sphere_transition(upper_bound_eps) == approx(1.0));
+    {
+      INFO("Extended");
+      const SphereTransition sphere_transition_extended{2., 4., true, true};
+      CHECK(sphere_transition_extended(lower_bound) == approx(0.0));
+      CHECK(sphere_transition_extended(lower_bound_eps) == approx(0.0));
+      const std::array<double, 3> inside_lower_bound{{1., 0., 0.}};
+      CHECK(sphere_transition_extended(inside_lower_bound) == 0.0);
+      CHECK(sphere_transition_extended(midpoint) == approx(0.5));
+      CHECK(sphere_transition_extended(upper_bound) == approx(1.));
+      CHECK(sphere_transition_extended(upper_bound_eps) == approx(1.));
+      const std::array<double, 3> outside_upper_bound{{5., 0., 0.}};
+      CHECK(sphere_transition_extended(outside_upper_bound) == 1.0);
+    }
   }
 }