From 59dbdfae51598acf9b72c5e5d7fb0e33ea5ea986 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 11 Mar 2024 20:49:25 +0100 Subject: [PATCH] Adjust interface --- src/CMakeLists.txt | 1 + src/utils/interpolation.hpp | 27 +++++++++++++++++---------- 2 files changed, 18 insertions(+), 10 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index de101a6f..05efed4d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -23,6 +23,7 @@ add_executable( utils/few_modes_ft.cpp utils/few_modes_ft.hpp utils/interpolation.hpp + utils/robust.hpp ) add_subdirectory(pgen) diff --git a/src/utils/interpolation.hpp b/src/utils/interpolation.hpp index 8cf17927..6185ffc9 100644 --- a/src/utils/interpolation.hpp +++ b/src/utils/interpolation.hpp @@ -1,4 +1,3 @@ - //======================================================================================== // AthenaPK - a performance portable block structured AMR astrophysical MHD code. // Copyright (c) 2024, Athena-Parthenon Collaboration. All rights reserved. @@ -23,23 +22,28 @@ #ifndef UTILS_INTERPOLATION_HPP_ #define UTILS_INTERPOLATION_HPP_ -// Spiner includes -#include - // Parthenon includes +#include "utils/error_checking.hpp" #include #include #include -// Phoebus includes -#include "phoebus_utils/grid_utils.hpp" -#include "phoebus_utils/robust.hpp" +#include "robust.hpp" namespace interpolation { using namespace parthenon::package::prelude; using parthenon::Coordinates_t; -using weights_t = Spiner::weights_t; + +// 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; + } +}; /// Base class for providing interpolation methods on uniformly spaced data. /// Constructor is provided with spacing, number of support points, and desired @@ -133,7 +137,10 @@ template KOKKOS_INLINE_FUNCTION void GetWeights(const Real x, const int nx, const Coordinates_t &coords, int &ix, weights_t &w) { - const Real min = Coordinates::GetXv(0, coords); + PARTHENON_DEBUG_REQUIRE( + typeid(Coordinates_t) == typeid(UniformCartesian), + "Interpolation routines currently only work for UniformCartesian"); + const Real min = coords.Xc(0); // assume uniform Cartesian const Real dx = coords.CellWidthFA(DIR); ix = std::min(std::max(0, static_cast(robust::ratio(x - min, dx))), nx - 2); const Real floor = min + ix * dx; @@ -235,4 +242,4 @@ KOKKOS_INLINE_FUNCTION Real Do(int b, const Real X1, const Real X2, const Real X // Convenience Namespace Alias namespace LCInterp = interpolation::Cent::Linear; -#endif // PHOEBUS_UTILS_PHOEBUS_INTERPOLATION_HPP_ +#endif // UTILS_INTERPOLATION_HPP_