diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 22a3bdac..6d106322 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -23,8 +23,6 @@ add_executable( tracers/tracers.hpp utils/few_modes_ft.cpp utils/few_modes_ft.hpp - utils/interpolation.hpp - utils/robust.hpp ) add_subdirectory(pgen) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 352361bc..1cd4dee1 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -23,17 +23,21 @@ #include #include -#include "../main.hpp" -#include "../utils/interpolation.hpp" +// Parthenon headers #include "basic_types.hpp" #include "interface/metadata.hpp" #include "kokkos_abstraction.hpp" #include "parthenon_array_generic.hpp" -#include "tracers.hpp" #include "utils/error_checking.hpp" +#include "utils/interpolation.hpp" + +// AthenaPK headers +#include "../main.hpp" +#include "tracers.hpp" namespace Tracers { using namespace parthenon::package::prelude; +namespace LCInterp = parthenon::interpolation::cent::linear; std::shared_ptr Initialize(ParameterInput *pin) { auto tracer_pkg = std::make_shared("tracers"); diff --git a/src/utils/interpolation.hpp b/src/utils/interpolation.hpp deleted file mode 100644 index 6185ffc9..00000000 --- a/src/utils/interpolation.hpp +++ /dev/null @@ -1,245 +0,0 @@ -//======================================================================================== -// AthenaPK - a performance portable block structured AMR astrophysical MHD code. -// Copyright (c) 2024, Athena-Parthenon Collaboration. All rights reserved. -// Licensed under the BSD 3-Clause License (the "LICENSE"). -//======================================================================================== -// 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_ - -// Parthenon includes -#include "utils/error_checking.hpp" -#include -#include -#include - -#include "robust.hpp" - -namespace interpolation { - -using namespace parthenon::package::prelude; -using parthenon::Coordinates_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 -/// shift. GetIndicesAndWeights then updates arrays of indices and weights for -/// calculating the interpolated data. These arrays are of size StencilSize(). -/// Data is forced to zero outside the boundaries. -class Interpolation { - public: - KOKKOS_FUNCTION - Interpolation(const int nSupport, const Real dx, const Real shift) - : nSupport_(nSupport), dx_(dx), shift_(shift), ishift_(std::round(shift)) {} - - KOKKOS_INLINE_FUNCTION - virtual void GetIndicesAndWeights(const int i, int *idx, Real *wgt) const {} - KOKKOS_INLINE_FUNCTION - virtual int StencilSize() const { return 0; } - - static constexpr int maxStencilSize = 2; - - protected: - const int nSupport_; - const Real dx_; - Real shift_; - int ishift_; -}; - -class PiecewiseConstant : public Interpolation { - public: - KOKKOS_FUNCTION - PiecewiseConstant(const int nSupport, const Real dx, const Real shift) - : Interpolation(nSupport, dx, shift) {} - - KOKKOS_INLINE_FUNCTION - void GetIndicesAndWeights(const int i, int *idx, Real *wgt) const override { - idx[0] = i + ishift_; - wgt[0] = 1.; - if (idx[0] < 0 || idx[0] >= nSupport_) { - idx[0] = 0; - wgt[0] = 0.; - } - } - - KOKKOS_INLINE_FUNCTION - int StencilSize() const override { return 1; } -}; - -class Linear : public Interpolation { - public: - KOKKOS_FUNCTION - Linear(const int nSupport, const Real dx, const Real shift) - : Interpolation(nSupport, dx, shift) { - PARTHENON_FAIL("Not written yet!"); - } - - KOKKOS_INLINE_FUNCTION - void GetIndicesAndWeights(const int i, int *idx, Real *wgt) const override { - idx[0] = std::floor(i + shift_); - idx[1] = idx[0] + 1; - - wgt[0] = wgt[1] = 1. - wgt[0]; - - for (int nsup = 0; nsup < 2; nsup++) { - if (idx[nsup] < 0 || idx[nsup] >= nSupport_) { - idx[nsup] = 0; - wgt[nsup] = 0.; - } - } - } - - KOKKOS_INLINE_FUNCTION - int StencilSize() const override { return 2; } -}; - -// 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 -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(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; - 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 -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(X1, p.GetDim(1), coords, ix[0], w[0]); - GetWeights(X2, p.GetDim(2), coords, ix[1], w[1]); - GetWeights(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 -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(X1, p.GetDim(1), coords, ix1, w1); - GetWeights(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 -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(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 -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 - -// Convenience Namespace Alias -namespace LCInterp = interpolation::Cent::Linear; - -#endif // UTILS_INTERPOLATION_HPP_ diff --git a/src/utils/robust.hpp b/src/utils/robust.hpp deleted file mode 100644 index 7488abea..00000000 --- a/src/utils/robust.hpp +++ /dev/null @@ -1,66 +0,0 @@ - -//======================================================================================== -// AthenaPK - a performance portable block structured AMR astrophysical MHD code. -// Copyright (c) 2024, Athena-Parthenon Collaboration. All rights reserved. -// Licensed under the BSD 3-Clause License (the "LICENSE"). -//======================================================================================== -// 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 "basic_types.hpp" -#include -#include - -namespace robust { -using parthenon::Real; - -template -KOKKOS_FORCEINLINE_FUNCTION constexpr auto LARGE() { - return 0.1 * std::numeric_limits::max(); -} -template -KOKKOS_FORCEINLINE_FUNCTION constexpr auto SMALL() { - return 10 * std::numeric_limits::min(); -} -template -KOKKOS_FORCEINLINE_FUNCTION constexpr auto EPS() { - return 10 * std::numeric_limits::epsilon(); -} - -template -KOKKOS_FORCEINLINE_FUNCTION auto make_positive(const T val) { - return std::max(val, EPS()); -} - -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 -KOKKOS_INLINE_FUNCTION int sgn(const T &val) { - return (T(0) <= val) - (val < T(0)); -} -template -KOKKOS_INLINE_FUNCTION auto ratio(const A &a, const B &b) { - const B sgn = b >= 0 ? 1 : -1; - return a / (b + sgn * SMALL()); -} -} // namespace robust - -#endif // UTILS_ROBUST_HPP_