Add bi- and trilinear interpolation (#1020)

* Add interpolation

* Make linter happy and remove dead code

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR1020]](https://github.com/parthenon-hpc-lab/parthenon/pull/1020) Add bi- and trilinear interpolation routines
 - [[PR 1026]](https://github.com/parthenon-hpc-lab/parthenon/pull/1026) Particle BCs without relocatable device code
 - [[PR 1037]](https://github.com/parthenon-hpc-lab/parthenon/pull/1037) Add SwarmPacks
 - [[PR 1068]](https://github.com/parthenon-hpc-lab/parthenon/pull/1068) Add ability to dump sparse pack contents as a string

benchmarks/burgers/recon.hpp

@@ -17,6 +17,7 @@
 #include <algorithm>
 
 #include "basic_types.hpp"
+#include <utils/robust.hpp>
 using parthenon::Real;
 
 namespace recon {
@@ -45,7 +46,7 @@ void WENO5Z(const Real q0, const Real q1, const Real q2, const Real q3, const Re
                                   {-1.0 / 6.0, 5.0 / 6.0, 1.0 / 3.0},
                                   {1.0 / 3.0, 5.0 / 6.0, -1.0 / 6.0}};
   constexpr Real w5gamma[3] = {0.1, 0.6, 0.3};
-  constexpr Real eps = 1e-100;
+  constexpr Real eps = parthenon::robust::EPS();
   constexpr Real thirteen_thirds = 13.0 / 3.0;
 
   Real a = q0 - 2 * q1 + q2;

src/CMakeLists.txt

@@ -253,6 +253,7 @@ add_library(parthenon
   utils/index_split.hpp
   utils/indexer.hpp
   utils/instrument.hpp
+  utils/interpolation.hpp
   utils/loop_utils.hpp
   utils/morton_number.hpp
   utils/mpi_types.hpp
@@ -261,6 +262,7 @@ add_library(parthenon
   utils/object_pool.hpp
   utils/partition_stl_containers.hpp
   utils/reductions.hpp
+  utils/robust.hpp
   utils/show_config.cpp
   utils/signal_handler.cpp
   utils/sort.hpp

src/utils/interpolation.hpp

@@ -0,0 +1,169 @@
+//========================================================================================
+// Parthenon performance portable AMR framework
+// Copyright(C) 2024 The Parthenon collaboration
+// Licensed under the 3-clause BSD License, see LICENSE file for details
+//========================================================================================
+// Interpolation copied/refactored from
+// https://github.com/lanl/phoebus and https://github.com/lanl/spiner
+//========================================================================================
+// © 2022. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
+// is operated by Triad National Security, LLC for the U.S.
+// Department of Energy/National Nuclear Security Administration. All
+// rights in the program are reserved by Triad National Security, LLC,
+// and the U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others
+// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works,
+// distribute copies to the public, perform publicly and display
+// publicly, and to permit others to do so.
+
+#ifndef UTILS_INTERPOLATION_HPP_
+#define UTILS_INTERPOLATION_HPP_
+
+#include <algorithm>
+
+// Parthenon includes
+#include <coordinates/coordinates.hpp>
+#include <kokkos_abstraction.hpp>
+#include <parthenon/package.hpp>
+#include <utils/error_checking.hpp>
+#include <utils/robust.hpp>
+
+namespace parthenon {
+namespace interpolation {
+
+// From https://github.com/lanl/spiner/blob/main/spiner/regular_grid_1d.hpp
+// a poor-man's std::pair
+struct weights_t {
+  Real first, second;
+  KOKKOS_INLINE_FUNCTION Real &operator[](const int i) {
+    assert(0 <= i && i <= 1);
+    return i == 0 ? first : second;
+  }
+};
+
+// TODO(JMM): Is this interpolation::Do syntax reasonable? An
+// alternative path would be a class called "LCInterp with all
+// static functions. Then it could have an `operator()` which would
+// be maybe nicer?
+// TODO(JMM): Merge this w/ what Ben has done.
+namespace cent {
+namespace linear {
+
+/*
+ * Get interpolation weights for linear interpolation
+ * PARAM[IN] - x - location to interpolate to
+ * PARAM[IN] - nx - number of points along this direction. Used for sanity checks.
+ * PARAM[IN] - coords - parthenon coords object
+ * PARAM[OUT] - ix - index of points to interpolate
+ * PARAM[OUT] - w - weights
+ */
+template <int DIR>
+KOKKOS_INLINE_FUNCTION void GetWeights(const Real x, const int nx,
+                                       const Coordinates_t &coords, int &ix,
+                                       weights_t &w) {
+  PARTHENON_DEBUG_REQUIRE(
+      typeid(Coordinates_t) == typeid(UniformCartesian),
+      "Interpolation routines currently only work for UniformCartesian");
+  const Real min = coords.Xc<DIR>(0); // assume uniform Cartesian
+  const Real dx = coords.CellWidthFA(DIR);
+  ix = std::min(std::max(0, static_cast<int>(robust::ratio(x - min, dx))), nx - 2);
+  const Real floor = min + ix * dx;
+  w[1] = robust::ratio(x - floor, dx);
+  w[0] = 1. - w[1];
+}
+
+/*
+ * Trilinear interpolation on a variable or meshblock pack
+ * PARAM[IN] - b - Meshblock index
+ * PARAM[IN] - X1, X2, X3 - Coordinate locations
+ * PARAM[IN] - p - Variable or MeshBlockPack
+ * PARAM[IN] - v - variable index
+ */
+template <typename Pack>
+KOKKOS_INLINE_FUNCTION Real Do3D(int b, const Real X1, const Real X2, const Real X3,
+                                 const Pack &p, int v) {
+  const auto &coords = p.GetCoords(b);
+  int ix[3];
+  weights_t w[3];
+  GetWeights<X1DIR>(X1, p.GetDim(1), coords, ix[0], w[0]);
+  GetWeights<X2DIR>(X2, p.GetDim(2), coords, ix[1], w[1]);
+  GetWeights<X3DIR>(X3, p.GetDim(3), coords, ix[2], w[2]);
+  return (w[2][0] * (w[1][0] * (w[0][0] * p(b, v, ix[2], ix[1], ix[0]) +
+                                w[0][1] * p(b, v, ix[2], ix[1], ix[0] + 1)) +
+                     w[1][1] * (w[0][0] * p(b, v, ix[2], ix[1] + 1, ix[0]) +
+                                w[0][1] * p(b, v, ix[2], ix[1] + 1, ix[0] + 1))) +
+          w[2][1] * (w[1][0] * (w[0][0] * p(b, v, ix[2] + 1, ix[1], ix[0]) +
+                                w[0][1] * p(b, v, ix[2] + 1, ix[1], ix[0] + 1)) +
+                     w[1][1] * (w[0][0] * p(b, v, ix[2] + 1, ix[1] + 1, ix[0]) +
+                                w[0][1] * p(b, v, ix[2] + 1, ix[1] + 1, ix[0] + 1))));
+}
+
+/*
+ * Bilinear interpolation on a variable or meshblock pack
+ * PARAM[IN] - b - Meshblock index
+ * PARAM[IN] - X1, X2 - Coordinate locations
+ * PARAM[IN] - p - Variable or MeshBlockPack
+ * PARAM[IN] - v - variable index
+ */
+template <typename Pack>
+KOKKOS_INLINE_FUNCTION Real Do2D(int b, const Real X1, const Real X2, const Pack &p,
+                                 int v) {
+  const auto &coords = p.GetCoords(b);
+  int ix1, ix2;
+  weights_t w1, w2;
+  GetWeights<X1DIR>(X1, p.GetDim(1), coords, ix1, w1);
+  GetWeights<X2DIR>(X2, p.GetDim(2), coords, ix2, w2);
+  return (w2[0] * (w1[0] * p(b, v, 0, ix2, ix1) + w1[1] * p(b, v, 0, ix2, ix1 + 1)) +
+          w2[1] *
+              (w1[0] * p(b, v, 0, ix2 + 1, ix1) + w1[1] * p(b, v, 0, ix2 + 1, ix1 + 1)));
+}
+
+/*
+ * Linear interpolation on a variable or meshblock pack
+ * PARAM[IN] - b - Meshblock index
+ * PARAM[IN] - X1 - Coordinate location
+ * PARAM[IN] - p - Variable or MeshBlockPack
+ * PARAM[IN] - v - variable index
+ */
+template <typename Pack>
+KOKKOS_INLINE_FUNCTION Real Do1D(int b, const Real X1, const Pack &p, int v) {
+  const auto &coords = p.GetCoords(b);
+  int ix;
+  weights_t w;
+  GetWeights<X1DIR>(X1, p.GetDim(1), coords, ix, w);
+  return w[0] * p(b, v, 0, 0, ix) + w[1] * p(b, v, 0, 0, ix + 1);
+}
+
+/*
+ * Trilinear or bilinear interpolation on a variable or meshblock pack
+ * PARAM[IN] - axisymmetric
+ * PARAM[IN] - b - Meshblock index
+ * PARAM[IN] - X1, X2, X3 - Coordinate locations
+ * PARAM[IN] - p - Variable or MeshBlockPack
+ * PARAM[IN] - v - variable index
+ */
+// JMM: I know this won't vectorize because of the switch, but it
+// probably won't anyway, since we're doing trilinear
+// interpolation, which will kill memory locality.  Doing it this
+// way means we can do trilinear vs bilinear which I think is a
+// sufficient win at minimum code bloat.
+template <typename Pack>
+KOKKOS_INLINE_FUNCTION Real Do(int b, const Real X1, const Real X2, const Real X3,
+                               const Pack &p, int v) {
+  if (p.GetDim(3) > 1) {
+    return Do3D(b, X1, X2, X3, p, v);
+  } else if (p.GetDim(2) > 1) {
+    return Do2D(b, X1, X2, p, v);
+  } else { // 1D
+    return Do1D(b, X1, p, v);
+  }
+}
+
+} // namespace linear
+} // namespace cent
+} // namespace interpolation
+} // namespace parthenon
+#endif // UTILS_INTERPOLATION_HPP_

src/utils/robust.hpp

@@ -0,0 +1,67 @@
+//========================================================================================
+// Parthenon performance portable AMR framework
+// Copyright(C) 2024 The Parthenon collaboration
+// Licensed under the 3-clause BSD License, see LICENSE file for details
+//========================================================================================
+// Copied from https://github.com/lanl/phoebus
+//========================================================================================
+// © 2022. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
+// is operated by Triad National Security, LLC for the U.S.
+// Department of Energy/National Nuclear Security Administration. All
+// rights in the program are reserved by Triad National Security, LLC,
+// and the U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others
+// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works,
+// distribute copies to the public, perform publicly and display
+// publicly, and to permit others to do so.
+
+#ifndef UTILS_ROBUST_HPP_
+#define UTILS_ROBUST_HPP_
+
+#include <algorithm>
+#include <limits>
+
+#include <config.hpp>
+#include <kokkos_abstraction.hpp>
+
+namespace parthenon {
+namespace robust {
+
+template <typename T = Real>
+KOKKOS_FORCEINLINE_FUNCTION constexpr auto LARGE() {
+  return 0.1 * std::numeric_limits<T>::max();
+}
+template <typename T = Real>
+KOKKOS_FORCEINLINE_FUNCTION constexpr auto SMALL() {
+  return 10 * std::numeric_limits<T>::min();
+}
+template <typename T = Real>
+KOKKOS_FORCEINLINE_FUNCTION constexpr auto EPS() {
+  return 10 * std::numeric_limits<T>::epsilon();
+}
+
+template <typename T>
+KOKKOS_FORCEINLINE_FUNCTION auto make_positive(const T val) {
+  return std::max(val, EPS<T>());
+}
+
+KOKKOS_FORCEINLINE_FUNCTION
+Real make_bounded(const Real val, const Real vmin, const Real vmax) {
+  return std::min(std::max(val, vmin + EPS()), vmax * (1.0 - EPS()));
+}
+
+template <typename T>
+KOKKOS_INLINE_FUNCTION int sgn(const T &val) {
+  return (T(0) <= val) - (val < T(0));
+}
+template <typename A, typename B>
+KOKKOS_INLINE_FUNCTION auto ratio(const A &a, const B &b) {
+  const B sgn = b >= 0 ? 1 : -1;
+  return a / (b + sgn * SMALL<B>());
+}
+} // namespace robust
+} // namespace parthenon
+#endif // UTILS_ROBUST_HPP_