Add alternate irregular interpolant

src/Evolution/DgSubcell/SetInterpolators.hpp

@@ -24,6 +24,7 @@
 #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
 #include "Evolution/DgSubcell/Tags/Reconstructor.hpp"
+#include "NumericalAlgorithms/Interpolation/ExtendedIrregularInterpolant.hpp"
 #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
 #include "NumericalAlgorithms/Spectral/Mesh.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"

src/NumericalAlgorithms/Interpolation/CMakeLists.txt

@@ -12,6 +12,7 @@ spectre_target_sources(
   BarycentricRationalSpanInterpolator.cpp
   CubicSpanInterpolator.cpp
   CubicSpline.cpp
+  ExtendedIrregularInterpolant.cpp
   IrregularInterpolant.cpp
   LinearLeastSquares.cpp
   LinearRegression.cpp
@@ -31,6 +32,7 @@ spectre_target_headers(
   BarycentricRationalSpanInterpolator.hpp
   CubicSpanInterpolator.hpp
   CubicSpline.hpp
+  ExtendedIrregularInterpolant.cpp
   IrregularInterpolant.hpp
   LagrangePolynomial.hpp
   LinearLeastSquares.hpp

src/NumericalAlgorithms/Interpolation/ExtendedIrregularInterpolant.cpp

@@ -0,0 +1,237 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "ExtendedIrregularInterpolant.hpp"
+
+#include <algorithm>
+#include <array>
+#include <iterator>
+#include <vector>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+// IWYU pragma: no_forward_declare Tensor
+
+namespace {
+
+// Just linear for now, can be extended to higher order...
+// Future optimization: it might be more efficient to use Blaze's sparse
+// matrices since the interpolation matrix is mostly zeros
+std::vector<double> fd_stencil(const DataVector& xi_source,
+                               const double xi_target) {
+  ASSERT(std::is_sorted(std::begin(xi_source), std::end(xi_source)),
+         "xi_source = " << xi_source);
+  auto xi_u =
+      std::upper_bound(std::begin(xi_source), std::end(xi_source), xi_target);
+  if (std::end(xi_source) == xi_u) {
+    std::advance(xi_u, -1);
+  }
+  if (std::begin(xi_source) == xi_u) {
+    std::advance(xi_u, 1);
+  }
+  const auto xi_l = std::prev(xi_u);
+  auto index = std::distance(std::begin(xi_source), xi_l);
+  std::vector<double> result(xi_source.size(), 0.0);
+  const auto result_l = std::next(std::begin(result), index);
+  const auto result_u = std::next(result_l, 1);
+  *result_l = (*xi_u - xi_target) / (*xi_u - *xi_l);
+  *result_u = (xi_target - *xi_l) / (*xi_u - *xi_l);
+  return result;
+}
+
+template <size_t Dim>
+Matrix interpolation_matrix(
+    const std::vector<size_t>& source_extents,
+    const tnsr::I<DataVector, Dim, Frame::ElementLogical>& source_points,
+    const tnsr::I<DataVector, Dim, Frame::ElementLogical>& target_points);
+
+template <>
+Matrix interpolation_matrix(
+    const std::vector<size_t>& source_extents,
+    const tnsr::I<DataVector, 1, Frame::ElementLogical>& source_points,
+    const tnsr::I<DataVector, 1, Frame::ElementLogical>& target_points) {
+  const auto number_of_source_points = source_extents[0];
+  const DataVector xi_source(
+      const_cast<DataVector&>(get<0>(source_points)).data(),
+      number_of_source_points);
+  const auto number_of_target_points = get<0>(target_points).size();
+  Matrix result(number_of_target_points, number_of_source_points);
+  for (size_t p = 0; p < number_of_target_points; ++p) {
+    const double xi_target = get<0>(target_points)[p];
+    const auto stencil = fd_stencil(xi_source, xi_target);
+    for (size_t i = 0; i < number_of_source_points; ++i) {
+      result(p, i) = stencil[i];
+    }
+  }
+  return result;
+}
+
+template <>
+Matrix interpolation_matrix(
+    const std::vector<size_t>& source_extents,
+    const tnsr::I<DataVector, 2, Frame::ElementLogical>& source_points,
+    const tnsr::I<DataVector, 2, Frame::ElementLogical>& target_points) {
+  const auto number_of_target_points = get<0>(target_points).size();
+  const auto number_of_source_points = get<0>(source_points).size();
+  Matrix result(number_of_target_points, number_of_source_points);
+
+  DataVector xi_source(const_cast<DataVector&>(get<0>(source_points)).data(),
+                       number_of_source_points);
+  DataVector eta_source;
+  if (source_extents[1] == source_extents[0]) {
+    eta_source.set_data_ref(&xi_source);
+  } else {
+    eta_source.destructive_resize(source_extents[1]);
+    for (size_t j = 0; j < source_extents[1]; ++j) {
+      eta_source[j] = get<1>(source_points)[j * source_extents[0]];
+    }
+  }
+    for (size_t p = 0; p < number_of_target_points; ++p) {
+      const double xi_target = get<0>(target_points)[p];
+      const double eta_target = get<1>(target_points)[p];
+      const auto xi_stencil = fd_stencil(xi_source, xi_target);
+      const auto eta_stencil = fd_stencil(eta_source, eta_target);
+      for (size_t j = 0, s = 0; j < source_extents[1]; ++j) {
+        for (size_t i = 0; i < source_extents[0]; ++i) {
+          result(p, s) = xi_stencil[i] * eta_stencil[j];
+          ++s;
+        }
+      }
+    }
+    return result;
+}
+
+template <>
+Matrix interpolation_matrix(
+    const std::vector<size_t>& source_extents,
+    const tnsr::I<DataVector, 3, Frame::ElementLogical>& source_points,
+    const tnsr::I<DataVector, 3, Frame::ElementLogical>& target_points) {
+  const auto number_of_target_points = get<0>(target_points).size();
+  const auto number_of_source_points = get<0>(source_points).size();
+
+  Matrix result(number_of_target_points, number_of_source_points);
+  DataVector xi_source(const_cast<DataVector&>(get<0>(source_points)).data(),
+                       number_of_source_points);
+  DataVector eta_source;
+  if (source_extents[1] == source_extents[0]) {
+    eta_source.set_data_ref(&xi_source);
+  } else {
+    eta_source.destructive_resize(source_extents[1]);
+    for (size_t j = 0; j < source_extents[1]; ++j) {
+      eta_source[j] = get<1>(source_points)[j * source_extents[0]];
+    }
+  }
+    DataVector zeta_source;
+    if (source_extents[2] == source_extents[0]) {
+      zeta_source.set_data_ref(&xi_source);
+    } else if (source_extents[2] == source_extents[1]) {
+      zeta_source.set_data_ref(&eta_source);
+    } else {
+      zeta_source.destructive_resize(source_extents[2]);
+      for (size_t k = 0; k < source_extents[2]; ++k) {
+        zeta_source[k] =
+            get<2>(source_points)[k * source_extents[0] * source_extents[1]];
+      }
+    }
+    for (size_t p = 0; p < number_of_target_points; ++p) {
+      const double xi_target = get<0>(target_points)[p];
+      const double eta_target = get<1>(target_points)[p];
+      const double zeta_target = get<2>(target_points)[p];
+      const auto xi_stencil = fd_stencil(xi_source, xi_target);
+      const auto eta_stencil = fd_stencil(eta_source, eta_target);
+      const auto zeta_stencil = fd_stencil(zeta_source, zeta_target);
+      for (size_t k = 0, s = 0; k < source_extents[2]; ++k) {
+        for (size_t j = 0; j < source_extents[1]; ++j) {
+          for (size_t i = 0; i < source_extents[0]; ++i) {
+            result(p, s) = xi_stencil[i] * eta_stencil[j] * zeta_stencil[k];
+            ++s;
+          }
+        }
+      }
+    }
+    return result;
+}
+}  // namespace
+
+namespace intrp {
+
+template <size_t Dim>
+ExtendedIrregular<Dim>::ExtendedIrregular() = default;
+
+template <size_t Dim>
+ExtendedIrregular<Dim>::ExtendedIrregular(
+    const std::vector<size_t> source_extents,
+    const tnsr::I<DataVector, Dim, Frame::ElementLogical>& source_points,
+    const tnsr::I<DataVector, Dim, Frame::ElementLogical>& target_points)
+    : interpolation_matrix_(
+          interpolation_matrix(source_extents, source_points, target_points)) {}
+
+template <size_t Dim>
+void ExtendedIrregular<Dim>::pup(PUP::er& p) {
+  p | interpolation_matrix_;
+}
+
+template <size_t Dim>
+void ExtendedIrregular<Dim>::interpolate(
+    const gsl::not_null<DataVector*> result, const DataVector& input) const {
+  const size_t m = interpolation_matrix_.rows();
+  const size_t k = interpolation_matrix_.columns();
+  ASSERT(k == input.size(),
+         "Number of points in 'input', "
+             << input.size()
+             << ",\n disagrees with the size of the source points, " << k
+             << ", that was passed into the constructor");
+  if (result->size() != m) {
+    result->destructive_resize(m);
+  }
+  dgemv_('n', m, k, 1.0, interpolation_matrix_.data(),
+         interpolation_matrix_.spacing(), input.data(), 1, 0.0, result->data(),
+         1);
+}
+
+template <size_t Dim>
+DataVector ExtendedIrregular<Dim>::interpolate(const DataVector& input) const {
+  DataVector result{input.size()};
+  interpolate(make_not_null(&result), input);
+  return result;
+}
+
+template <size_t Dim>
+void ExtendedIrregular<Dim>::interpolate(
+    const gsl::not_null<gsl::span<double>*> result,
+    const gsl::span<const double>& input) const {
+  const size_t m = interpolation_matrix_.rows();
+  const size_t k = interpolation_matrix_.columns();
+  ASSERT(input.size() % k == 0,
+         "Number of points in 'input', "
+             << input.size()
+             << ",\n must be a multiple of the source grid points, " << k
+             << ", that was passed into the constructor");
+  const size_t number_of_components = input.size() / k;
+  ASSERT(result->size() == number_of_components * m,
+         "The result must be of size " << number_of_components * m
+                                       << " but got " << result->size());
+  dgemm_<true>('N', 'N', interpolation_matrix_.rows(), number_of_components,
+               interpolation_matrix_.columns(), 1.0,
+               interpolation_matrix_.data(), interpolation_matrix_.rows(),
+               input.data(), interpolation_matrix_.columns(), 0.0,
+               result->data(), interpolation_matrix_.rows());
+}
+
+template <size_t Dim>
+bool operator!=(const ExtendedIrregular<Dim>& lhs,
+                const ExtendedIrregular<Dim>& rhs) {
+  return not(lhs == rhs);
+}
+
+template class ExtendedIrregular<1>;
+template class ExtendedIrregular<2>;
+template class ExtendedIrregular<3>;
+template bool operator!=(const ExtendedIrregular<1>& lhs,
+                         const ExtendedIrregular<1>& rhs);
+template bool operator!=(const ExtendedIrregular<2>& lhs,
+                         const ExtendedIrregular<2>& rhs);
+template bool operator!=(const ExtendedIrregular<3>& lhs,
+                         const ExtendedIrregular<3>& rhs);
+
+}  // namespace intrp

src/NumericalAlgorithms/Interpolation/ExtendedIrregularInterpolant.hpp

@@ -0,0 +1,119 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Matrix.hpp"
+#include "DataStructures/Tensor/TypeAliases.hpp"
+#include "DataStructures/Variables.hpp"
+#include "Utilities/Blas.hpp"
+#include "Utilities/ErrorHandling/Assert.hpp"
+#include "Utilities/Gsl.hpp"
+
+/// \cond
+template <size_t Dim>
+class Mesh;
+namespace PUP {
+class er;
+}  // namespace PUP
+// IWYU pragma: no_forward_declare Variables
+/// \endcond
+
+namespace intrp {
+
+/// \ingroup NumericalAlgorithmsGroup
+/// \brief Interpolates a `Variables` onto an arbitrary set of points.
+///
+/// \details This is an alternate version of IrregularInterpolant, taking
+/// the extents and source points as inputs instead of the source mesh.
+/// The purpose of this class is for particular cases in FD where
+/// we were getting problematic extrapolation (Sphere domain) for ghost point
+/// reconstruction in the case where the target points lied outside
+/// of the nearest neighbor.
+template <size_t Dim>
+class ExtendedIrregular {
+ public:
+  ExtendedIrregular(
+      const std::vector<size_t> source_extents,
+      const tnsr::I<DataVector, Dim, Frame::ElementLogical>& source_points,
+      const tnsr::I<DataVector, Dim, Frame::ElementLogical>& target_points);
+  ExtendedIrregular();
+
+  /// Serialization for Charm++
+  // NOLINTNEXTLINE(google-runtime-references)
+  void pup(PUP::er& p);
+
+  /// @{
+  /// Performs the interpolation on a `Variables` with grid points corresponding
+  /// to the source grid points as specified in the constructor.
+  /// The result is a `Variables` whose internal `DataVector` goes over the
+  /// list of target_points that were specified in the constructor.
+  /// \note for the void function, `result` will be resized to the proper size.
+  template <typename TagsList>
+  void interpolate(gsl::not_null<Variables<TagsList>*> result,
+                   const Variables<TagsList>& vars) const;
+  template <typename TagsList>
+  Variables<TagsList> interpolate(const Variables<TagsList>& vars) const;
+  /// @}
+
+  /// @{
+  /// \brief Interpolate a DataVector onto the target points.
+  ///
+  /// \note When interpolating multiple tensors, the Variables interface is more
+  /// efficient. However, this DataVector interface is useful for applications
+  /// where only some components of a Tensor or Variables need to be
+  /// interpolated.
+  void interpolate(gsl::not_null<DataVector*> result,
+                   const DataVector& input) const;
+  DataVector interpolate(const DataVector& input) const;
+  /// @}
+
+  /// \brief Interpolate multiple variables on the grid to the target points.
+  void interpolate(gsl::not_null<gsl::span<double>*> result,
+                   const gsl::span<const double>& input) const;
+
+ private:
+  friend bool operator==(const ExtendedIrregular& lhs,
+                         const ExtendedIrregular& rhs) {
+    return lhs.interpolation_matrix_ == rhs.interpolation_matrix_;
+  }
+  Matrix interpolation_matrix_;
+};
+
+template <size_t Dim>
+template <typename TagsList>
+void ExtendedIrregular<Dim>::interpolate(
+    const gsl::not_null<Variables<TagsList>*> result,
+    const Variables<TagsList>& vars) const {
+  if (UNLIKELY(result->number_of_grid_points() !=
+               interpolation_matrix_.rows())) {
+    *result = Variables<TagsList>(interpolation_matrix_.rows(), 0.);
+  }
+  ASSERT(interpolation_matrix_.columns() == vars.number_of_grid_points(),
+         "Number of grid points in source 'vars', "
+             << vars.number_of_grid_points()
+             << ",\n disagrees with the size of the source_mesh, "
+             << interpolation_matrix_.columns()
+             << ", that was passed into the constructor");
+  gsl::span<double> result_span{result->data(), result->size()};
+  const gsl::span<const double> vars_span{vars.data(), vars.size()};
+  interpolate(make_not_null(&result_span), vars_span);
+}
+
+template <size_t Dim>
+template <typename TagsList>
+Variables<TagsList> ExtendedIrregular<Dim>::interpolate(
+    const Variables<TagsList>& vars) const {
+  Variables<TagsList> result{interpolation_matrix_.rows()};
+  interpolate(make_not_null(&result), vars);
+  return result;
+}
+
+template <size_t Dim>
+bool operator!=(const ExtendedIrregular<Dim>& lhs,
+                const ExtendedIrregular<Dim>& rhs);
+
+}  // namespace intrp