Exporter: add mode to extrapolate from nearest element

Uses the Lagrange basis in the nearest element to extrapolate into the
excision region. This can be unstable when we want to extrapolate
far into the excision because we're extrapolating over many logical element sizes.
However, it can be useful for extrapolating only very slightly into the
excision, e.g. if the excision in the initial data is the apparent horizon
and the evolution needs a little space around it to find and track the
horizon.

src/IO/Exporter/CMakeLists.txt

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

src/IO/Exporter/Exporter.cpp

@@ -28,6 +28,59 @@
 #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 +382,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 +499,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 +574,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 +611,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))

src/IO/Exporter/Exporter.hpp

@@ -15,6 +15,14 @@
 #include <variant>
 #include <vector>
 
+/// \cond
+namespace Options {
+class Option;
+template <typename T>
+struct create_from_yaml;
+}  // namespace Options
+/// \endcond
+
 /// Functions that are intended to be used by external programs, e.g. to
 /// interpolate data in volume files to target points.
 namespace spectre::Exporter {
@@ -32,6 +40,46 @@ struct ObservationStep {
   int value;
 };
 
+enum class ExcisionExtrapolationMode {
+  /*!
+   * \brief No extrapolation into the excisions is performed
+   */
+  NoExtrapolation,
+  /*!
+   * \brief Extrapolate into excisions from the nearest element
+   *
+   * Uses the Lagrange basis in the nearest element to extrapolate into the
+   * excision region. This can be unstable when the element is too small (e.g.
+   * when it is h-refined a lot) and/or has log or inverse radial map. In those
+   * cases the logical coordinate that we're extrapolating to can be many
+   * logical element sizes away. Also, this relies on inverting the neighboring
+   * block's map outside the block, which is not guaranteed to work.
+   *
+   * This mode can be useful for extrapolating only very slightly into the
+   * excision, e.g. if the excision in the initial data is the apparent horizon
+   * and the evolution needs a little space around it to find and track the
+   * horizon.
+   */
+  NearestElement,
+  /*!
+   * \brief Extrapolate into excisions from radial anchor points
+   *
+   * Uses the strategy from \cite Etienne2008re (adjusted for distorted
+   * excisions): we choose uniformly spaced radial anchor points spaced as
+   * $\Delta r = 0.3 r_\mathrm{AH}$ in the grid frame (where the excision is
+   * spherical), then map the anchor points to the distorted frame (where we
+   * have the target point) and do a 7th order polynomial extrapolation into the
+   * excision region.
+   *
+   * This mode is the most stable and accurate to fill the entire excision for
+   * moving puncture codes.
+   */
+  RadialAnchors,
+};
+
+std::ostream& operator<<(std::ostream& os,
+                         const ExcisionExtrapolationMode& mode);
+
 /*!
  * \brief Interpolate data in volume files to target points
  *
@@ -47,15 +95,12 @@ struct ObservationStep {
  * "Lapse", "Shift_x", "Shift_y", "Shift_z", "SpatialMetric_xx", etc.
  * Look into the H5 file to see what components are available.
  * \param target_points The points to interpolate to, in inertial coordinates.
- * \param extrapolate_into_excisions Enables extrapolation into excision regions
- * of the domain (default is `false`). This can be useful to fill the excision
- * region with (constraint-violating but smooth) data so it can be imported into
- * moving puncture codes. Specifically, we implement the strategy used in
- * \cite Etienne2008re adjusted for distorted excisions: we choose uniformly
- * spaced radial anchor points spaced as $\Delta r = 0.3 r_\mathrm{AH}$ in the
- * grid frame (where the excision is spherical), then map the anchor points to
- * the distorted frame (where we have the target point) and do a 7th order
- * polynomial extrapolation into the excision region.
+ * \param excision_extrapolation_mode Enables extrapolation into excision
+ * regions of the domain (default is `NoExtrapolation`). This can be useful to
+ * fill the excision region with (constraint-violating but smooth) data so it
+ * can be imported into moving puncture codes, or to extrapolate a little bit
+ * into the excision to help the evolution find and track the apparent horizon.
+ * See `ExcisionExtrapolationMode` for the available modes.
  * \param num_threads The number of threads to use if OpenMP is linked in. If
  * not specified, OpenMP will determine the number of threads automatically.
  * It's also possible to set the number of threads using the environment
@@ -73,7 +118,22 @@ 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,
-    bool extrapolate_into_excisions = false,
+    ExcisionExtrapolationMode excision_extrapolation_mode =
+        ExcisionExtrapolationMode::NoExtrapolation,
     std::optional<size_t> num_threads = std::nullopt);
 
 }  // namespace spectre::Exporter
+
+template <>
+struct Options::create_from_yaml<spectre::Exporter::ExcisionExtrapolationMode> {
+  template <typename Metavariables>
+  static spectre::Exporter::ExcisionExtrapolationMode create(
+      const Options::Option& options) {
+    return create<void>(options);
+  }
+};
+
+template <>
+spectre::Exporter::ExcisionExtrapolationMode
+Options::create_from_yaml<spectre::Exporter::ExcisionExtrapolationMode>::create<
+    void>(const Options::Option& options);

src/IO/Exporter/InterpolateToPoints.hpp

@@ -37,7 +37,8 @@ tuples::tagged_tuple_from_typelist<Tags> interpolate_to_points(
     const std::string& subfile_name,
     const std::variant<ObservationId, ObservationStep>& observation,
     const tnsr::I<DataType, Dim>& target_points,
-    bool extrapolate_into_excisions = false,
+    ExcisionExtrapolationMode excision_extrapolation_mode =
+        ExcisionExtrapolationMode::NoExtrapolation,
     std::optional<size_t> num_threads = std::nullopt) {
   // Convert target_points to an array of vectors
   // Possible performance optimization: avoid copying the data by allowing
@@ -65,7 +66,7 @@ tuples::tagged_tuple_from_typelist<Tags> interpolate_to_points(
   // Interpolate!
   const auto interpolated_data = interpolate_to_points(
       volume_files_or_glob, subfile_name, observation, tensor_components,
-      target_points_array, extrapolate_into_excisions, num_threads);
+      target_points_array, excision_extrapolation_mode, num_threads);
   // Convert the interpolated data to a tagged_tuple
   tuples::tagged_tuple_from_typelist<Tags> result{};
   size_t component_index = 0;

src/IO/Exporter/Python/Bindings.cpp

@@ -14,32 +14,41 @@ namespace py = pybind11;
 
 PYBIND11_MODULE(_Pybindings, m) {  // NOLINT
   enable_segfault_handler();
+  py::enum_<spectre::Exporter::ExcisionExtrapolationMode>(
+      m, "ExcisionExtrapolationMode")
+      .value("NoExtrapolation",
+             spectre::Exporter::ExcisionExtrapolationMode::NoExtrapolation)
+      .value("NearestElement",
+             spectre::Exporter::ExcisionExtrapolationMode::NearestElement)
+      .value("RadialAnchors",
+             spectre::Exporter::ExcisionExtrapolationMode::RadialAnchors);
   m.def(
       "interpolate_to_points",
       [](const std::variant<std::vector<std::string>, std::string>&
              volume_files_or_glob,
          const std::string& subfile_name, const size_t observation_id,
          const std::vector<std::string>& tensor_components,
          std::vector<std::vector<double>> target_points,
-         bool extrapolate_into_excisions,
+         spectre::Exporter::ExcisionExtrapolationMode
+             excision_extrapolation_mode,
          const std::optional<size_t>& num_threads) {
         const size_t dim = target_points.size();
         const spectre::Exporter::ObservationId obs_id{observation_id};
         if (dim == 1) {
           return spectre::Exporter::interpolate_to_points(
               volume_files_or_glob, subfile_name, obs_id, tensor_components,
               make_array<std::vector<double>, 1>(std::move(target_points)),
-              extrapolate_into_excisions, num_threads);
+              excision_extrapolation_mode, num_threads);
         } else if (dim == 2) {
           return spectre::Exporter::interpolate_to_points(
               volume_files_or_glob, subfile_name, obs_id, tensor_components,
               make_array<std::vector<double>, 2>(std::move(target_points)),
-              extrapolate_into_excisions, num_threads);
+              excision_extrapolation_mode, num_threads);
         } else if (dim == 3) {
           return spectre::Exporter::interpolate_to_points(
               volume_files_or_glob, subfile_name, obs_id, tensor_components,
               make_array<std::vector<double>, 3>(std::move(target_points)),
-              extrapolate_into_excisions, num_threads);
+              excision_extrapolation_mode, num_threads);
         } else {
           ERROR("Invalid dimension of target points: "
                 << dim
@@ -50,6 +59,8 @@ PYBIND11_MODULE(_Pybindings, m) {  // NOLINT
       },
       py::arg("volume_files_or_glob"), py::arg("subfile_name"),
       py::arg("observation_id"), py::arg("tensor_components"),
-      py::arg("target_points"), py::arg("extrapolate_into_excisions") = false,
+      py::arg("target_points"),
+      py::arg("excision_extrapolation_mode") =
+          spectre::Exporter::ExcisionExtrapolationMode::NoExtrapolation,
       py::arg("num_threads") = std::nullopt);
 }

src/IO/Exporter/Python/InterpolateToPoints.py

@@ -8,7 +8,7 @@
 import rich
 
 import spectre.IO.H5 as spectre_h5
-from spectre.IO.Exporter import interpolate_to_points
+from spectre.IO.Exporter import ExcisionExtrapolationMode, interpolate_to_points
 from spectre.Visualization.OpenVolfiles import (
     open_volfiles_command,
     parse_points,
@@ -43,12 +43,20 @@
 )
 @click.option(
     "--extrapolate-into-excisions",
-    is_flag=True,
+    "excision_extrapolation_mode",
+    type=click.Choice(ExcisionExtrapolationMode.__members__),
+    default=ExcisionExtrapolationMode.NoExtrapolation.name,
+    show_default=True,
+    callback=lambda ctx, param, value: ExcisionExtrapolationMode.__members__[
+        value
+    ],
     help=(
-        "Enables extrapolation into excision regions of the domain. "
-        "This can be useful to fill the excision region with "
-        "(constraint-violating but smooth) data so it can be imported into "
-        "moving puncture codes."
+        "Enables extrapolation into excision regions of the domain. This can be"
+        " useful to fill the excision region with (constraint-violating but"
+        " smooth) data so it can be imported into moving puncture codes (set to"
+        " 'AnchorPoints'), or to extrapolate a little bit into the excision to"
+        " help the evolution find and track the apparent horizon (set to"
+        " 'NearestElement')."
     ),
 )
 @click.option(

src/PointwiseFunctions/InitialDataUtilities/NumericData.cpp

@@ -9,18 +9,19 @@
 #include <string>
 #include <utility>
 
-NumericData::NumericData(std::string file_glob, std::string subgroup,
-                         int observation_step, bool extrapolate_into_excisions)
+NumericData::NumericData(
+    std::string file_glob, std::string subgroup, int observation_step,
+    spectre::Exporter::ExcisionExtrapolationMode excision_extrapolation_mode)
     : file_glob_(std::move(file_glob)),
       subgroup_(std::move(subgroup)),
       observation_step_(observation_step),
-      extrapolate_into_excisions_(extrapolate_into_excisions) {}
+      excision_extrapolation_mode_(excision_extrapolation_mode) {}
 
 bool operator==(const NumericData& lhs, const NumericData& rhs) {
   return lhs.file_glob() == rhs.file_glob() and
          lhs.subgroup() == rhs.subgroup() and
          lhs.observation_step() == rhs.observation_step() and
-         lhs.extrapolate_into_excisions() == rhs.extrapolate_into_excisions();
+         lhs.excision_extrapolation_mode() == rhs.excision_extrapolation_mode();
 }
 
 bool operator!=(const NumericData& lhs, const NumericData& rhs) {
@@ -31,7 +32,7 @@ void NumericData::pup(PUP::er& p) {
   p | file_glob_;
   p | subgroup_;
   p | observation_step_;
-  p | extrapolate_into_excisions_;
+  p | excision_extrapolation_mode_;
 }
 
 namespace elliptic::analytic_data {