Extrapolate from nearest element

src/Domain/BlockLogicalCoordinates.cpp

@@ -21,7 +21,8 @@ template <size_t Dim, typename Fr>
 std::optional<tnsr::I<double, Dim, ::Frame::BlockLogical>>
 block_logical_coordinates_single_point(
     const tnsr::I<double, Dim, Fr>& input_point, const Block<Dim>& block,
-    const double time, const domain::FunctionsOfTimeMap& functions_of_time) {
+    const double time, const domain::FunctionsOfTimeMap& functions_of_time,
+    const bool allow_extrapolation) {
   std::optional<tnsr::I<double, Dim, ::Frame::BlockLogical>> logical_point{};
   if (block.is_time_dependent()) {
     if constexpr (std::is_same_v<Fr, ::Frame::Inertial>) {
@@ -125,14 +126,57 @@ block_logical_coordinates_single_point(
       logical_point->get(d) = -1.0;
       continue;
     }
-    if (abs(logical_point->get(d)) > 1.0) {
+    if (abs(logical_point->get(d)) > 1.0 and not allow_extrapolation) {
       return std::nullopt;
     }
   }
 
   return logical_point;
 }
 
+template <size_t Dim, typename Fr>
+BlockLogicalCoords<Dim> block_logical_coordinates_in_excision(
+    const tnsr::I<double, Dim, Fr>& input_point,
+    const ExcisionSphere<Dim>& excision_sphere,
+    const std::vector<Block<Dim>>& blocks, const double time,
+    const domain::FunctionsOfTimeMap& functions_of_time) {
+  // Comment by NV (Apr 2024): This function can be made more robust by first
+  // checking if the point is inside the excision sphere at all, but the
+  // excision sphere doesn't currently have this information. It needs the
+  // grid-to-inertial map including the deformation of the excision sphere to
+  // determine this.
+  for (const auto& [block_id, direction] :
+       excision_sphere.abutting_directions()) {
+    auto x_logical = block_logical_coordinates_single_point(
+        input_point, blocks[block_id], time, functions_of_time, true);
+    if (not x_logical.has_value()) {
+      continue;
+    }
+    // Discard block if the point has angular logical coordinates outside the
+    // range [-1, 1]
+    for (size_t d = 0; d < Dim; ++d) {
+      if (d != direction.dimension() and std::abs(x_logical->get(d)) > 1.) {
+        x_logical = std::nullopt;
+        break;
+      }
+    }
+    if (not x_logical.has_value()) {
+      continue;
+    }
+    // Discard block if the point is radially inside the block or on the other
+    // side of the excision sphere
+    const double radial_distance_this_block =
+        x_logical->get(direction.dimension()) * direction.sign();
+    if (radial_distance_this_block < 1.) {
+      continue;
+    }
+    // The checks above should leave only 1 valid block, so return that
+    return make_id_pair(domain::BlockId(block_id),
+                        std::move(x_logical.value()));
+  }
+  return std::nullopt;
+}
+
 template <size_t Dim, typename Fr>
 std::vector<BlockLogicalCoords<Dim>> block_logical_coordinates(
     const Domain<Dim>& domain, const tnsr::I<DataVector, Dim, Fr>& x,
@@ -174,11 +218,18 @@ std::vector<BlockLogicalCoords<Dim>> block_logical_coordinates(
   block_logical_coordinates_single_point(                                      \
       const tnsr::I<double, DIM(data), FRAME(data)>& input_point,              \
       const Block<DIM(data)>& block, const double time,                        \
-      const domain::FunctionsOfTimeMap& functions_of_time);                    \
+      const domain::FunctionsOfTimeMap& functions_of_time,                     \
+      bool allow_extrapolation);                                               \
   template std::vector<BlockLogicalCoords<DIM(data)>>                          \
   block_logical_coordinates(                                                   \
       const Domain<DIM(data)>& domain,                                         \
       const tnsr::I<DataVector, DIM(data), FRAME(data)>& x, const double time, \
+      const domain::FunctionsOfTimeMap& functions_of_time);                    \
+  template BlockLogicalCoords<DIM(data)>                                       \
+  block_logical_coordinates_in_excision(                                       \
+      const tnsr::I<double, DIM(data), FRAME(data)>& input_point,              \
+      const ExcisionSphere<DIM(data)>& excision_sphere,                        \
+      const std::vector<Block<DIM(data)>>& blocks, const double time,          \
       const domain::FunctionsOfTimeMap& functions_of_time);
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3),

src/Domain/BlockLogicalCoordinates.hpp

@@ -21,6 +21,8 @@ template <size_t VolumeDim>
 class Domain;
 template <size_t VolumeDim>
 class Block;
+template <size_t VolumeDim>
+class ExcisionSphere;
 /// \endcond
 
 template <size_t Dim>
@@ -69,10 +71,32 @@ auto block_logical_coordinates(
     const domain::FunctionsOfTimeMap& functions_of_time = {})
     -> std::vector<BlockLogicalCoords<Dim>>;
 
+/// \par Extrapolation
+/// If `allow_extrapolation` is `true` (default is `false`), then the function
+/// can return block-logical coordinates outside of [-1, 1]. This is useful for
+/// extrapolating into excision regions.
 template <size_t Dim, typename Fr>
 std::optional<tnsr::I<double, Dim, ::Frame::BlockLogical>>
 block_logical_coordinates_single_point(
     const tnsr::I<double, Dim, Fr>& input_point, const Block<Dim>& block,
     double time = std::numeric_limits<double>::signaling_NaN(),
-    const domain::FunctionsOfTimeMap& functions_of_time = {});
+    const domain::FunctionsOfTimeMap& functions_of_time = {},
+    bool allow_extrapolation = false);
 /// @}
+
+/*!
+ * \brief Finds the block-logical coordinates of a point in an excision sphere.
+ *
+ * The point will be reported to belong to the nearest block that abuts the
+ * excision sphere and will have a logical coordinate outside the range [-1, 1]
+ * in the radial direction. This is useful for extrapolating into excision
+ * regions using the Lagrange polynomial basis of the nearest element. If the
+ * point is not within the excision sphere, the return value is std::nullopt.
+ */
+template <size_t Dim, typename Frame>
+BlockLogicalCoords<Dim> block_logical_coordinates_in_excision(
+    const tnsr::I<double, Dim, Frame>& input_point,
+    const ExcisionSphere<Dim>& excision_sphere,
+    const std::vector<Block<Dim>>& blocks,
+    double time = std::numeric_limits<double>::signaling_NaN(),
+    const domain::FunctionsOfTimeMap& functions_of_time = {});

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

@@ -65,10 +65,12 @@ std::optional<double> SphereTransition::original_radius_over_radius(
   if ((original_radius + eps_) >= r_min_ and
       (original_radius - eps_) <= r_max_) {
     return std::optional<double>{original_radius / mag};
-  } else if ((a_ > 0.0 and mag > r_max_) or (a_ < 0.0 and mag < r_min_)) {
+  } else if ((a_ > 0.0 and original_radius > r_max_) or
+             (a_ < 0.0 and original_radius < r_min_)) {
     // a_ being positive is a sentinel for reverse (see constructor)
     return {1.0 + radial_distortion / mag};
-  } else if ((a_ < 0.0 and mag > r_max_) or (a_ > 0.0 and mag < r_min_)) {
+  } else if ((a_ < 0.0 and original_radius > r_max_) or
+             (a_ > 0.0 and original_radius < r_min_)) {
     return {1.0};
   } else {
     return std::nullopt;

src/Domain/ElementLogicalCoordinates.cpp

@@ -27,36 +27,29 @@ element_logical_coordinates(
     const ElementId<Dim>& element_id) {
   tnsr::I<double, Dim, Frame::ElementLogical> x_element_logical{};
   for (size_t d = 0; d < Dim; ++d) {
-    // Check if the point is outside the element
-    const double up = element_id.segment_id(d).endpoint(Side::Upper);
-    const double lo = element_id.segment_id(d).endpoint(Side::Lower);
-    if (x_block_logical.get(d) < lo or x_block_logical.get(d) > up) {
+    // Check if the point is inside the element in this dimension.
+    // The following conditions are valid:
+    // - The block-logical coordinates are within the element's bounds.
+    // - The block-logical coordinate is outside the block and we're in the
+    //   outermost element in that direction. This means we want to extrapolate
+    //   out of the block.
+    const auto& segment_id = element_id.segment_id(d);
+    const double up = segment_id.endpoint(Side::Upper);
+    const double lo = segment_id.endpoint(Side::Lower);
+    if ((x_block_logical.get(d) >= lo and x_block_logical.get(d) <= up) or
+        (x_block_logical.get(d) < -1. and segment_id.index() == 0) or
+        (x_block_logical.get(d) > 1. and
+         segment_id.index() == two_to_the(segment_id.refinement_level()) - 1)) {
+      // Map to element logical coords
+      x_element_logical.get(d) =
+          (2.0 * x_block_logical.get(d) - up - lo) / (up - lo);
+    } else {
       return std::nullopt;
     }
-    // Map to element logical coords
-    x_element_logical.get(d) =
-        (2.0 * x_block_logical.get(d) - up - lo) / (up - lo);
   }
   return x_element_logical;
 }
 
-namespace {
-// The segments bounds are binary fractions (i.e. the numerator is an
-// integer and the denominator is a power of 2) so these floating point
-// comparisons should be safe from roundoff problems
-// Need to return true if on upper face of block
-// Otherwise return true if point is within the segment or on the lower bound
-bool segment_contains(const double x_block_logical,
-                      const double lower_bound_block_logical,
-                      const double upper_bound_block_logical) {
-  if (UNLIKELY(x_block_logical == upper_bound_block_logical)) {
-    return (upper_bound_block_logical == 1.0);
-  }
-  return (x_block_logical >= lower_bound_block_logical and
-          x_block_logical < upper_bound_block_logical);
-}
-}  // namespace
-
 template <size_t Dim>
 std::unordered_map<ElementId<Dim>, ElementLogicalCoordHolder<Dim>>
 element_logical_coordinates(
@@ -90,13 +83,15 @@ element_logical_coordinates(
         if (not x_elem.has_value()) {
           continue;
         }
-        // Disambiguate points on shared element boundaries
+        // Disambiguate points on shared element boundaries. To do this,
+        // associate boundary points with the element with the lower index.
         bool is_contained = true;
         for (size_t d = 0; d < Dim; ++d) {
-          const double up = element_id.segment_id(d).endpoint(Side::Upper);
-          const double lo = element_id.segment_id(d).endpoint(Side::Lower);
-          const double x_block_log = x_block_logical.get(d);
-          if (not segment_contains(x_block_log, lo, up)) {
+          const auto& segment_id = element_id.segment_id(d);
+          const double up = segment_id.endpoint(Side::Upper);
+          if (x_block_logical.get(d) == up and
+              segment_id.index() !=
+                  two_to_the(segment_id.refinement_level()) - 1) {
             is_contained = false;
             break;
           }

src/Domain/ElementLogicalCoordinates.hpp

@@ -32,6 +32,11 @@ class IdPair;
  * have element logical coordinates of -1 or 1. See the other function overload
  * for handling multiple points and disambiguating points on shared element
  * boundaries.
+ *
+ * Points with block-logical coordinates outside [-1, 1] are considered to
+ * belong to the outermost element of the block in that direction, and will
+ * therefore return an element logical coordinate outside of [-1, 1] as well.
+ * This is useful for extrapolation.
  */
 template <size_t Dim>
 std::optional<tnsr::I<double, Dim, Frame::ElementLogical>>

src/IO/Exporter/CMakeLists.txt

@@ -34,6 +34,7 @@ target_link_libraries(
   FunctionsOfTime
   H5
   Interpolation
+  Options
   Serialization
   Utilities
   )

src/IO/Exporter/Exporter.cpp

@@ -21,13 +21,68 @@
 #include "IO/H5/VolumeData.hpp"
 #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
 #include "NumericalAlgorithms/Interpolation/PolynomialInterpolation.hpp"
+#include "Options/Context.hpp"
+#include "Options/Options.hpp"
 #include "Utilities/FileSystem.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/GetOutput.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/Overloader.hpp"
 #include "Utilities/Serialization/Serialize.hpp"
 
+namespace spectre::Exporter {
+
+std::ostream& operator<<(std::ostream& os,
+                         const ExcisionExtrapolationMode& mode) {
+  switch (mode) {
+    case ExcisionExtrapolationMode::NoExtrapolation:
+      os << "NoExtrapolation";
+      break;
+    case ExcisionExtrapolationMode::NearestElement:
+      os << "NearestElement";
+      break;
+    case ExcisionExtrapolationMode::RadialAnchors:
+      os << "RadialAnchors";
+      break;
+    default:
+      os << "Unknown";
+      break;
+  }
+  return os;
+}
+
+}  // namespace spectre::Exporter
+
+template <>
+spectre::Exporter::ExcisionExtrapolationMode
+Options::create_from_yaml<spectre::Exporter::ExcisionExtrapolationMode>::create<
+    void>(const Options::Option& options) {
+  const auto type_read = options.parse_as<std::string>();
+  for (const auto value :
+       {spectre::Exporter::ExcisionExtrapolationMode::NoExtrapolation,
+        spectre::Exporter::ExcisionExtrapolationMode::NearestElement,
+        spectre::Exporter::ExcisionExtrapolationMode::RadialAnchors}) {
+    if (type_read == get_output(value)) {
+      return value;
+    }
+  }
+  PARSE_ERROR(
+      options.context(),
+      "Failed to convert \""
+          << type_read
+          << "\" to spectre::Exporter::ExcisionExtrapolationMode.\nMust be one "
+             "of "
+          << get_output(
+                 spectre::Exporter::ExcisionExtrapolationMode::NoExtrapolation)
+          << ", "
+          << get_output(
+                 spectre::Exporter::ExcisionExtrapolationMode::NearestElement)
+          << ", or "
+          << get_output(
+                 spectre::Exporter::ExcisionExtrapolationMode::RadialAnchors)
+          << ".");
+}
+
 // Ignore OpenMP pragmas when OpenMP is not enabled
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Wunknown-pragmas"
@@ -329,7 +384,7 @@ std::vector<std::vector<double>> interpolate_to_points(
     const std::variant<ObservationId, ObservationStep>& observation,
     const std::vector<std::string>& tensor_components,
     const std::array<std::vector<double>, Dim>& target_points,
-    const bool extrapolate_into_excisions,
+    const ExcisionExtrapolationMode excision_extrapolation_mode,
     const std::optional<size_t> num_threads) {
   domain::creators::register_derived_with_charm();
   domain::creators::time_dependence::register_derived_with_charm();
@@ -446,19 +501,37 @@ std::vector<std::vector<double>> interpolate_to_points(
         }
       }  // for blocks
       if (block_logical_coords[s].has_value() or
-          not extrapolate_into_excisions) {
+          excision_extrapolation_mode ==
+              ExcisionExtrapolationMode::NoExtrapolation) {
         continue;
       }
-      // The point wasn't found in any block. Check if it's in an excision and
-      // set up extrapolation if requested.
+      // The point wasn't found in any block. Extrapolate if requested.
       for (const auto& [name, excision_sphere] : domain.excision_spheres()) {
-        if (add_extrapolation_anchors(
-                make_not_null(&extra_block_logical_coords),
-                make_not_null(&extra_extrapolation_info), excision_sphere,
-                domain, target_point, time, functions_of_time,
-                extrapolation_spacing)) {
-          extra_extrapolation_info.back().target_index = s;
-          break;
+        if (excision_extrapolation_mode ==
+            ExcisionExtrapolationMode::NearestElement) {
+          // Get the block-logical coordinates for the nearest element (the
+          // coordinates will be outside [-1, 1], so Lagrange extrapolation will
+          // happen)
+          block_logical_coords[s] = block_logical_coordinates_in_excision(
+              target_point, excision_sphere, domain.blocks(), time,
+              functions_of_time);
+          if (block_logical_coords[s].has_value()) {
+            break;
+          }
+        } else if (excision_extrapolation_mode ==
+                   ExcisionExtrapolationMode::RadialAnchors) {
+          // Set up extrapolation anchors for this point
+          if (add_extrapolation_anchors(
+                  make_not_null(&extra_block_logical_coords),
+                  make_not_null(&extra_extrapolation_info), excision_sphere,
+                  domain, target_point, time, functions_of_time,
+                  extrapolation_spacing)) {
+            extra_extrapolation_info.back().target_index = s;
+            break;
+          }
+        } else {
+          ERROR("Invalid excision extrapolation mode: "
+                << excision_extrapolation_mode);
         }
       }
     }  // omp for target points
@@ -503,7 +576,7 @@ std::vector<std::vector<double>> interpolate_to_points(
     }
   }
 
-  if (extrapolate_into_excisions) {
+  if (excision_extrapolation_mode == ExcisionExtrapolationMode::RadialAnchors) {
     // Extrapolate into excisions from the anchor points
 #pragma omp parallel for num_threads(resolved_num_threads)
     for (const auto& extrapolation : extrapolation_info) {
@@ -540,7 +613,7 @@ std::vector<std::vector<double>> interpolate_to_points(
       const std::variant<ObservationId, ObservationStep>& observation,        \
       const std::vector<std::string>& tensor_components,                      \
       const std::array<std::vector<double>, DIM(data)>& target_points,        \
-      bool extrapolate_into_excisions,                                        \
+      ExcisionExtrapolationMode excision_extrapolation_mode,                  \
       const std::optional<size_t> num_threads);
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))