From 72fa4565f210774fff67a2bb19cf9bffe2c23d2a Mon Sep 17 00:00:00 2001 From: Jooheon Yoo Date: Wed, 26 Jun 2024 12:53:23 -0400 Subject: [PATCH] Add alternate irregular interpolant --- src/Evolution/DgSubcell/SetInterpolators.hpp | 1 + .../Interpolation/CMakeLists.txt | 2 + .../ExtendedIrregularInterpolant.cpp | 237 ++++++++++++++++++ .../ExtendedIrregularInterpolant.hpp | 119 +++++++++ 4 files changed, 359 insertions(+) create mode 100644 src/NumericalAlgorithms/Interpolation/ExtendedIrregularInterpolant.cpp create mode 100644 src/NumericalAlgorithms/Interpolation/ExtendedIrregularInterpolant.hpp diff --git a/src/Evolution/DgSubcell/SetInterpolators.hpp b/src/Evolution/DgSubcell/SetInterpolators.hpp index 7852ae34f15e5..7e7ca42698abb 100644 --- a/src/Evolution/DgSubcell/SetInterpolators.hpp +++ b/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" diff --git a/src/NumericalAlgorithms/Interpolation/CMakeLists.txt b/src/NumericalAlgorithms/Interpolation/CMakeLists.txt index a9f9bc4bcd6e5..5a3fb864e97aa 100644 --- a/src/NumericalAlgorithms/Interpolation/CMakeLists.txt +++ b/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 diff --git a/src/NumericalAlgorithms/Interpolation/ExtendedIrregularInterpolant.cpp b/src/NumericalAlgorithms/Interpolation/ExtendedIrregularInterpolant.cpp new file mode 100644 index 0000000000000..3c4523148bd4e --- /dev/null +++ b/src/NumericalAlgorithms/Interpolation/ExtendedIrregularInterpolant.cpp @@ -0,0 +1,237 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "ExtendedIrregularInterpolant.hpp" + +#include +#include +#include +#include + +#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 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 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 +Matrix interpolation_matrix( + const std::vector& source_extents, + const tnsr::I& source_points, + const tnsr::I& target_points); + +template <> +Matrix interpolation_matrix( + const std::vector& source_extents, + const tnsr::I& source_points, + const tnsr::I& target_points) { + const auto number_of_source_points = source_extents[0]; + const DataVector xi_source( + const_cast(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& source_extents, + const tnsr::I& source_points, + const tnsr::I& 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(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& source_extents, + const tnsr::I& source_points, + const tnsr::I& 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(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 +ExtendedIrregular::ExtendedIrregular() = default; + +template +ExtendedIrregular::ExtendedIrregular( + const std::vector source_extents, + const tnsr::I& source_points, + const tnsr::I& target_points) + : interpolation_matrix_( + interpolation_matrix(source_extents, source_points, target_points)) {} + +template +void ExtendedIrregular::pup(PUP::er& p) { + p | interpolation_matrix_; +} + +template +void ExtendedIrregular::interpolate( + const gsl::not_null 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 +DataVector ExtendedIrregular::interpolate(const DataVector& input) const { + DataVector result{input.size()}; + interpolate(make_not_null(&result), input); + return result; +} + +template +void ExtendedIrregular::interpolate( + const gsl::not_null*> result, + const gsl::span& 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_('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 +bool operator!=(const ExtendedIrregular& lhs, + const ExtendedIrregular& 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 diff --git a/src/NumericalAlgorithms/Interpolation/ExtendedIrregularInterpolant.hpp b/src/NumericalAlgorithms/Interpolation/ExtendedIrregularInterpolant.hpp new file mode 100644 index 0000000000000..0645290804698 --- /dev/null +++ b/src/NumericalAlgorithms/Interpolation/ExtendedIrregularInterpolant.hpp @@ -0,0 +1,119 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#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 +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 +class ExtendedIrregular { + public: + ExtendedIrregular( + const std::vector source_extents, + const tnsr::I& source_points, + const tnsr::I& 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 + void interpolate(gsl::not_null*> result, + const Variables& vars) const; + template + Variables interpolate(const Variables& 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 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*> result, + const gsl::span& input) const; + + private: + friend bool operator==(const ExtendedIrregular& lhs, + const ExtendedIrregular& rhs) { + return lhs.interpolation_matrix_ == rhs.interpolation_matrix_; + } + Matrix interpolation_matrix_; +}; + +template +template +void ExtendedIrregular::interpolate( + const gsl::not_null*> result, + const Variables& vars) const { + if (UNLIKELY(result->number_of_grid_points() != + interpolation_matrix_.rows())) { + *result = Variables(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 result_span{result->data(), result->size()}; + const gsl::span vars_span{vars.data(), vars.size()}; + interpolate(make_not_null(&result_span), vars_span); +} + +template +template +Variables ExtendedIrregular::interpolate( + const Variables& vars) const { + Variables result{interpolation_matrix_.rows()}; + interpolate(make_not_null(&result), vars); + return result; +} + +template +bool operator!=(const ExtendedIrregular& lhs, + const ExtendedIrregular& rhs); + +} // namespace intrp