Exporter: allow extrapolation into excisions

docs/References.bib

@@ -843,6 +843,21 @@ @article{Dumbser2017okk
     year = "2018"
 }
 
+@article{Etienne2008re,
+  author        = {Etienne, Zachariah B. and Liu, Yuk Tung and Shapiro, Stuart
+                   L. and Baumgarte, Thomas W.},
+  title         = {General relativistic simulations of black-hole-neutron-star
+                   mergers: Effects of black-hole spin},
+  eprint        = {0812.2245},
+  archiveprefix = {arXiv},
+  primaryclass  = {astro-ph},
+  doi           = {10.1103/PhysRevD.79.044024},
+  journal       = {Phys. Rev. D},
+  volume        = {79},
+  pages         = {044024},
+  year          = {2009}
+}
+
 @article{Etienne2010ui,
   author        = "Etienne, Zachariah B. and Liu, Yuk Tung and Shapiro, Stuart
                   L.",

src/DataStructures/Tensor/EagerMath/CartesianToSpherical.cpp

@@ -85,7 +85,7 @@ tnsr::I<DataType, Dim, CoordsFrame> spherical_to_cartesian(
                                        Frame::Spherical<FRAME(data)>>& x);
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector), (1, 2, 3),
-                        (Frame::Grid, Frame::Inertial))
+                        (Frame::Grid, Frame::Distorted, Frame::Inertial))
 
 #undef DTYPE
 #undef DIM

src/IO/Exporter/Exporter.cpp

@@ -8,6 +8,7 @@
 #include <omp.h>
 #endif  // _OPENMP
 
+#include "DataStructures/Tensor/EagerMath/CartesianToSpherical.hpp"
 #include "Domain/BlockLogicalCoordinates.hpp"
 #include "Domain/Creators/RegisterDerivedWithCharm.hpp"
 #include "Domain/Creators/TimeDependence/RegisterDerivedWithCharm.hpp"
@@ -19,6 +20,7 @@
 #include "IO/H5/TensorData.hpp"
 #include "IO/H5/VolumeData.hpp"
 #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
+#include "NumericalAlgorithms/Interpolation/PolynomialInterpolation.hpp"
 #include "Utilities/FileSystem.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/GetOutput.hpp"
@@ -123,7 +125,176 @@ void interpolate_to_points(
         }
       }
     }  // omp for
-  }    // omp parallel
+  }  // omp parallel
+}
+
+// Data structure for extrapolation of tensor components into excisions
+template <size_t NumExtrapolationAnchors>
+struct ExtrapolationInfo {
+  // Index of the target point in the result array (y target)
+  size_t target_index;
+  // Index of the first anchor point in the source array (y data)
+  size_t source_index;
+  // Coordinate of the target point (x target)
+  double target_point;
+  // Coordinates of the anchor points (x data)
+  std::array<double, NumExtrapolationAnchors> anchors;
+};
+
+// Construct anchor points for extrapolation into excisions.
+//
+// Anchor points are placed radially in the grid frame around the excision
+// sphere, which is spherical in the grid frame. Then the anchor points are
+// added to the `block_logical_coords` so data is interpolated to them. Also an
+// entry is added to `extrapolation_info` so the interpolated data on the anchor
+// points can be extrapolated to the target point in the excision.
+//
+// This function returns `true` if the target point is inside the excision and
+// the anchor points were added. Otherwise, returns `false` and the next
+// excision should be tried.
+//
+// Alternatives and notes:
+// - We could extrapolate directly from the nearest element, using the Lagrange
+//   polynomial basis that the element already has. However, this is unstable
+//   when the element is too small (possibly h-refined) and/or has log spacing.
+//   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.
+// - We could construct anchor points in different frames, e.g. in the inertial
+//   frame. These choices probably don't make a big difference.
+// - We could do more performance optimizations for the case where the excision
+//   sphere is spherical. E.g. in that case the radial anchor points for all
+//   points are the same.
+template <size_t Dim, size_t NumExtrapolationAnchors>
+bool add_extrapolation_anchors(
+    const gsl::not_null<std::vector<BlockLogicalCoords<Dim>>*>
+        block_logical_coords,
+    const gsl::not_null<
+        std::vector<ExtrapolationInfo<NumExtrapolationAnchors>>*>
+        extrapolation_info,
+    const ExcisionSphere<Dim>& excision_sphere, const Domain<Dim>& domain,
+    const tnsr::I<double, Dim, Frame::Inertial>& target_point,
+    const double time, const domain::FunctionsOfTimeMap& functions_of_time,
+    const double extrapolation_spacing) {
+  // Get spherical coordinates around the excision in distorted frame.
+  // Note that the excision sphere doesn't have a distorted map, but its
+  // grid-to-inertial map is just the neighboring block's
+  // distorted-to-inertial map, so we can use either of those.
+  auto x_distorted =
+      [&excision_sphere, &target_point, &time,
+       &functions_of_time]() -> tnsr::I<double, Dim, Frame::Distorted> {
+    if (excision_sphere.is_time_dependent()) {
+      const auto x_grid =
+          excision_sphere.moving_mesh_grid_to_inertial_map().inverse(
+              target_point, time, functions_of_time);
+      ASSERT(x_grid.has_value(),
+             "Failed to invert grid-to-inertial map of excision sphere at "
+                 << excision_sphere.center() << " for point " << target_point
+                 << ".");
+      tnsr::I<double, Dim, Frame::Distorted> result{};
+      for (size_t d = 0; d < Dim; ++d) {
+        result.get(d) = x_grid->get(d);
+      }
+      return result;
+    } else {
+      tnsr::I<double, Dim, Frame::Distorted> result{};
+      for (size_t d = 0; d < Dim; ++d) {
+        result.get(d) = target_point.get(d);
+      }
+      return result;
+    }
+  }();
+  for (size_t d = 0; d < Dim; ++d) {
+    x_distorted.get(d) -= excision_sphere.center().get(d);
+  }
+  auto x_spherical_distorted = cartesian_to_spherical(x_distorted);
+  // Construct anchor points in grid frame
+  tnsr::I<DataVector, Dim, Frame::Spherical<Frame::Grid>>
+      anchors_grid_spherical{NumExtrapolationAnchors};
+  for (size_t i = 0; i < NumExtrapolationAnchors; ++i) {
+    get<0>(anchors_grid_spherical)[i] =
+        1. + static_cast<double>(i) * extrapolation_spacing;
+  }
+  get<0>(anchors_grid_spherical) *= excision_sphere.radius();
+  // The grid-distorted map preserves angles
+  if constexpr (Dim > 1) {
+    get<1>(anchors_grid_spherical) = get<1>(x_spherical_distorted);
+  }
+  if constexpr (Dim > 2) {
+    get<2>(anchors_grid_spherical) = get<2>(x_spherical_distorted);
+  }
+  auto anchors_grid = spherical_to_cartesian(anchors_grid_spherical);
+  for (size_t d = 0; d < Dim; ++d) {
+    anchors_grid.get(d) += excision_sphere.center().get(d);
+  }
+  // Map anchor points to block logical coordinates. These are needed for
+  // interpolation.
+  auto anchors_block_logical =
+      block_logical_coordinates(domain, anchors_grid, time, functions_of_time);
+  const size_t block_id = anchors_block_logical[0]->id.get_index();
+  const auto& block = domain.blocks()[block_id];
+  // Map anchor points to the distorted frame. We will extrapolate in
+  // the distorted frame because we have the target point in the distorted
+  // frame. The target point in the grid frame is undefined because there's
+  // no grid-distorted map in the excision sphere. We could extrapolate in
+  // the inertial frame, but that's just an unnecessary transformation.
+  auto anchors_distorted =
+      [&block, &anchors_grid, &time,
+       &functions_of_time]() -> tnsr::I<DataVector, Dim, Frame::Distorted> {
+    if (block.has_distorted_frame()) {
+      return block.moving_mesh_grid_to_distorted_map()(anchors_grid, time,
+                                                       functions_of_time);
+    } else {
+      tnsr::I<DataVector, Dim, Frame::Distorted> result{};
+      for (size_t d = 0; d < Dim; ++d) {
+        result.get(d) = anchors_grid.get(d);
+      }
+      return result;
+    }
+  }();
+  for (size_t d = 0; d < Dim; ++d) {
+    anchors_distorted.get(d) -= excision_sphere.center().get(d);
+  }
+  // Compute magnitude(anchors_distorted) and store in std::array
+  std::array<double, NumExtrapolationAnchors> radial_anchors_distorted_frame{};
+  for (size_t i = 0; i < NumExtrapolationAnchors; ++i) {
+    auto& radial_anchor_distorted_frame =
+        gsl::at(radial_anchors_distorted_frame, i);
+    radial_anchor_distorted_frame = square(get<0>(anchors_distorted)[i]);
+    if constexpr (Dim > 1) {
+      radial_anchor_distorted_frame += square(get<1>(anchors_distorted)[i]);
+    }
+    if constexpr (Dim > 2) {
+      radial_anchor_distorted_frame += square(get<2>(anchors_distorted)[i]);
+    }
+    radial_anchor_distorted_frame = sqrt(radial_anchor_distorted_frame);
+  }
+  // Return false if the target point is not inside the excision. It would be
+  // nice to do this earlier to avoid unnecessary work. Here we use the first
+  // anchor point in the distorted frame, which is at the excision boundary in
+  // the angular direction of the target point.
+  const double excision_radius_distorted = radial_anchors_distorted_frame[0];
+  if (get<0>(x_spherical_distorted) > excision_radius_distorted) {
+    return false;
+  }
+  // Add the anchor points to the the interpolation target points
+  for (size_t i = 0; i < NumExtrapolationAnchors; ++i) {
+    ASSERT(anchors_block_logical[i].has_value(),
+           "Extrapolation anchor point is not in any block. This should "
+           "not happen.");
+    block_logical_coords->push_back(std::move(anchors_block_logical[i]));
+  }
+  // Add the anchor points to the extrapolation info
+  extrapolation_info->push_back(
+      {// Target index is the point we're extrapolating to. It will be set
+       // outside this function.
+       std::numeric_limits<size_t>::max(),
+       // Source index will be adjusted outside this function as well.
+       extrapolation_info->size() * NumExtrapolationAnchors,
+       // Target point and anchors are the radii in the distorted frame
+       get<0>(x_spherical_distorted),
+       std::move(radial_anchors_distorted_frame)});
+  return true;
 }
 
 // Determines the selected observation ID in the volume data file, given either
@@ -158,6 +329,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 std::optional<size_t> num_threads) {
   domain::creators::register_derived_with_charm();
   domain::creators::time_dependence::register_derived_with_charm();
@@ -235,10 +407,21 @@ std::vector<std::vector<double>> interpolate_to_points(
   // through the domain. This is the most expensive part of the function, so we
   // parallelize the loop.
   std::vector<BlockLogicalCoords<Dim>> block_logical_coords(num_target_points);
+  // We also set up the extrapolation into excisions here. Anchor points are
+  // added to the `block_logical_coords` and additional information is collected
+  // in `extrapolation_info` for later extrapolation.
+  constexpr size_t num_extrapolation_anchors = 8;
+  const double extrapolation_spacing = 0.3;
+  std::vector<ExtrapolationInfo<num_extrapolation_anchors>>
+      extrapolation_info{};
 #pragma omp parallel num_threads(resolved_num_threads)
   {
+    // Set up thread-local variables
     tnsr::I<double, Dim, Frame::Inertial> target_point{};
-#pragma omp for
+    std::vector<BlockLogicalCoords<Dim>> extra_block_logical_coords{};
+    std::vector<ExtrapolationInfo<num_extrapolation_anchors>>
+        extra_extrapolation_info{};
+#pragma omp for nowait
     for (size_t s = 0; s < num_target_points; ++s) {
       for (size_t d = 0; d < Dim; ++d) {
         target_point.get(d) = gsl::at(target_points, d)[s];
@@ -252,17 +435,49 @@ std::vector<std::vector<double>> interpolate_to_points(
           break;
         }
       }  // for blocks
-    }    // omp for target points
-  }      // omp parallel
+      if (block_logical_coords[s].has_value() or
+          not extrapolate_into_excisions) {
+        continue;
+      }
+      // The point wasn't found in any block. Check if it's in an excision and
+      // set up extrapolation 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;
+        }
+      }
+    }  // omp for target points
+#pragma omp critical
+    {
+      // Append the extra block logical coordinates and extrapolation info from
+      // this thread to the global vectors. Also set the source index to the
+      // index in `block_logical_coords` where we're going to insert the new
+      // coordinates.
+      for (auto& info : extra_extrapolation_info) {
+        info.source_index += block_logical_coords.size();
+      }
+      block_logical_coords.insert(block_logical_coords.end(),
+                                  extra_block_logical_coords.begin(),
+                                  extra_block_logical_coords.end());
+      extrapolation_info.insert(extrapolation_info.end(),
+                                extra_extrapolation_info.begin(),
+                                extra_extrapolation_info.end());
+    }  // omp critical
+  }  // omp parallel
 
   // Allocate memory for result
   std::vector<std::vector<double>> result{};
   result.reserve(tensor_components.size());
   for (size_t i = 0; i < tensor_components.size(); ++i) {
-    result.emplace_back(num_target_points,
+    result.emplace_back(block_logical_coords.size(),
                         std::numeric_limits<double>::signaling_NaN());
   }
-  std::vector<bool> filled_data(num_target_points, false);
+  std::vector<bool> filled_data(block_logical_coords.size(), false);
 
   // Process all volume files in serial, because loading data with H5 must be
   // done in serial anyway. Instead, the loop over elements within each file is
@@ -277,6 +492,29 @@ std::vector<std::vector<double>> interpolate_to_points(
       break;
     }
   }
+
+  if (extrapolate_into_excisions) {
+    // Extrapolate into excisions from the anchor points
+#pragma omp parallel for num_threads(resolved_num_threads)
+    for (const auto& extrapolation : extrapolation_info) {
+      double extrapolation_error = 0.;
+      for (size_t i = 0; i < tensor_components.size(); ++i) {
+        intrp::polynomial_interpolation<num_extrapolation_anchors - 1>(
+            make_not_null(&result[i][extrapolation.target_index]),
+            make_not_null(&extrapolation_error), extrapolation.target_point,
+            // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
+            gsl::make_span(result[i].data() + extrapolation.source_index,
+                           num_extrapolation_anchors),
+            gsl::make_span(extrapolation.anchors.data(),
+                           num_extrapolation_anchors));
+      }
+    }
+    // Clear the anchor points from the result
+    for (size_t i = 0; i < tensor_components.size(); ++i) {
+      result[i].resize(num_target_points);
+    }
+  }
+
   return result;
 }
 
@@ -292,6 +530,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,                                        \
       const std::optional<size_t> num_threads);
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))

src/IO/Exporter/Exporter.hpp

@@ -45,6 +45,15 @@ 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 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
@@ -62,6 +71,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,
+    bool extrapolate_into_excisions = false,
     std::optional<size_t> num_threads = std::nullopt);
 
 }  // namespace spectre::Exporter