coordinates.hh

Go to the documentation of this file.

#pragma once
 
// standard library
#include <cmath>
#include <cstddef>
#include <cstdio>
#include <cstdlib>
#include <limits>
#include <stdexcept>
#include <tuple>
#include <type_traits>
#include <vector>
 
// DiFfRG
#include <DiFfRG/common/kokkos.hh>
#include <DiFfRG/common/utils.hh>
 
namespace DiFfRG
{
  // A concept for what is a coordinate class
  template <typename T>
  concept is_coordinates = requires(T t) {
    typename T::ctype;
    T::dim;
    t.size();
    t.forward(device::array<size_t, T::dim>{});
    { t.from_linear_index(size_t{}) } -> std::same_as<device::array<size_t, T::dim>>;
  };
 
  namespace internal
  {
    template <typename T, typename = void> struct coord_periodic : std::false_type {
    };
    template <typename T>
    struct coord_periodic<T, std::void_t<decltype(T::periodic)>> : std::bool_constant<T::periodic> {
    };
  } // namespace internal
 
  template <typename T> inline constexpr bool is_periodic_coordinate_v = internal::coord_periodic<T>::value;
 
  namespace internal
  {
    template <typename C, size_t i, typename = void> struct axis_periodic : std::false_type {
    };
    template <typename C, size_t i>
    struct axis_periodic<C, i, std::void_t<typename C::template coordinate_type<i>>>
        : std::bool_constant<is_periodic_coordinate_v<typename C::template coordinate_type<i>>> {
    };
  } // namespace internal
 
  template <typename C, size_t i> inline constexpr bool is_periodic_axis_v = internal::axis_periodic<C, i>::value;
 
  template <typename... Coordinates> class CoordinatePackND
  {
    // Assert that all coordinates have the same ctype
    static_assert((std::is_same<typename Coordinates::ctype,
                                typename device::tuple_element<0, device::tuple<Coordinates...>>::type::ctype>::value &&
                   ...),
                  "CoordinatePackND: all coordinate systems must use the same ctype");
    static_assert(sizeof...(Coordinates) > 0, "CoordinatePackND requires at least one coordinate system");
    // Assert that all coordinates have the dim 1
    static_assert(((Coordinates::dim == 1) && ...), "CoordinatePackND requires all coordinates to have dim 1");
 
  public:
    using ctype = typename device::tuple_element<0, device::tuple<Coordinates...>>::type::ctype;
    static constexpr size_t dim = sizeof...(Coordinates);
 
    CoordinatePackND(Coordinates... coordinates) : coordinates(coordinates...) {}
 
    template <typename... I>
      requires(std::is_convertible_v<std::decay_t<I>, size_t> && ...)
    KOKKOS_FORCEINLINE_FUNCTION device::array<ctype, dim> forward(I &&...i) const
    {
      static_assert(sizeof...(I) == sizeof...(Coordinates));
      return forward_impl(device::make_integer_sequence<int, sizeof...(I)>(), device::forward<I>(i)...);
    }
    template <typename IT>
    KOKKOS_FORCEINLINE_FUNCTION device::array<ctype, dim> forward(const device::array<IT, dim> &coords) const
    {
      return device::apply([&](const auto &...iargs) { return forward(iargs...); }, coords);
    }
 
    template <typename... I, int... Is>
    device::array<ctype, dim> KOKKOS_FORCEINLINE_FUNCTION forward_impl(device::integer_sequence<int, Is...>,
                                                                       I &&...i) const
    {
      return {{device::get<Is>(coordinates).forward(device::get<Is>(device::tie(i...)))...}};
    }
 
    template <typename... I> KOKKOS_FORCEINLINE_FUNCTION device::array<ctype, sizeof...(I)> backward(I &&...i) const
    {
      static_assert(sizeof...(I) == sizeof...(Coordinates));
      return backward_impl(device::make_integer_sequence<int, sizeof...(I)>(), device::forward<I>(i)...);
    }
 
    template <typename... I, int... Is>
    device::array<ctype, sizeof...(I)> KOKKOS_FORCEINLINE_FUNCTION backward_impl(device::integer_sequence<int, Is...>,
                                                                                 I &&...i) const
    {
      return {{device::get<Is>(coordinates).backward(device::get<Is>(device::tie(i...)))...}};
    }
 
    template <size_t i> const auto &get_coordinates() const { return device::get<i>(coordinates); }
 
    template <size_t i>
    using coordinate_type = typename device::tuple_element<i, device::tuple<Coordinates...>>::type;
 
    size_t KOKKOS_FORCEINLINE_FUNCTION size() const
    {
      // multiply all sizes
      size_t size = 1;
      constexpr_for<0, sizeof...(Coordinates), 1>([&](auto i) { size *= device::get<i>(coordinates).size(); });
      return size;
    }
 
    device::array<size_t, sizeof...(Coordinates)> KOKKOS_FORCEINLINE_FUNCTION sizes() const
    {
      device::array<size_t, sizeof...(Coordinates)> sizes;
      constexpr_for<0, sizeof...(Coordinates), 1>([&](auto i) { sizes[i] = device::get<i>(coordinates).size(); });
      return sizes;
    }
 
    device::array<size_t, sizeof...(Coordinates)> KOKKOS_INLINE_FUNCTION from_linear_index(size_t s) const
    {
      device::array<size_t, sizeof...(Coordinates)> idx;
      // calculate the index for each coordinate system
      constexpr_for<0, sizeof...(Coordinates), 1>([&](auto i) {
        idx[sizeof...(Coordinates) - 1 - i] = s % device::get<sizeof...(Coordinates) - 1 - i>(coordinates).size();
        s = s / device::get<sizeof...(Coordinates) - 1 - i>(coordinates).size();
      });
 
      return idx;
    }
 
    template <typename... Coordinates2>
    friend bool operator==(const CoordinatePackND<Coordinates...> &r, const CoordinatePackND<Coordinates2...> &l)
    {
      if constexpr (sizeof...(Coordinates) != sizeof...(Coordinates2)) return false; // Different number of coordinates
      if constexpr ((!std::is_same_v<Coordinates, Coordinates2> && ...))
        return false; // Different types of coordinates
      else if constexpr (sizeof...(Coordinates) == 0)
        return true; // Empty coordinate pack, always equal
      else {
        // Check if all coordinates are equal
        bool equal = true;
        constexpr_for<0, sizeof...(Coordinates), 1>(
            [&](auto i) { equal &= (device::get<i>(r.coordinates) == device::get<i>(l.coordinates)); });
        return equal;
      }
    }
 
    std::string to_string() const
    {
      std::string ret = "CoordinatePackND(";
      constexpr_for<0, dim, 1>([&](auto i) {
        ret += device::get<i>(coordinates).to_string();
        if (i + 1 < dim) ret += ", ";
      });
      ret += ")";
      return ret;
    }
 
  protected:
    const device::tuple<Coordinates...> coordinates;
  };
 
  template <typename Base> class SubCoordinates : public Base
  {
  public:
    using ctype = typename Base::ctype;
    static constexpr size_t dim = Base::dim;
 
    SubCoordinates(const Base &base, size_t offset, size_t size) : Base(base), m_offset(offset), m_size(size)
    {
      if (size == 0) throw std::runtime_error("SubCoordinates: size must be > 0");
      if (offset + size > base.size()) throw std::runtime_error("SubCoordinates: offset + size must be <= base.size()");
    }
 
    device::array<size_t, dim> KOKKOS_FORCEINLINE_FUNCTION sizes() const
    {
      static_assert(dim == 1, "SubCoordinates is a linear index window, which is a per-axis box only for dim == 1.");
      return device::array<size_t, dim>{{m_size}};
    }
 
    size_t KOKKOS_FORCEINLINE_FUNCTION size() const { return m_size; }
 
    device::array<size_t, dim> KOKKOS_INLINE_FUNCTION from_linear_index(size_t s) const
    {
      return Base::from_linear_index(m_offset + s);
    }
 
    template <typename... I>
      requires(std::is_convertible_v<std::decay_t<I>, size_t> && ...)
    KOKKOS_FORCEINLINE_FUNCTION device::array<ctype, dim> forward(I &&...i) const
    {
      static_assert(sizeof...(I) == dim);
      return forward(device::array<size_t, dim>{{static_cast<size_t>(i)...}});
    }
 
    template <typename IT> KOKKOS_INLINE_FUNCTION device::array<ctype, dim> forward(device::array<IT, dim> i) const
    {
      return Base::forward(i);
    }
 
    // backward() is deliberately not overridden: it is the inverse of forward(), which is now
    // Base's unchanged, so Base::backward is already the right answer and is inherited as is.
 
    std::string to_string() const
    {
      return "SubCoordinates(" + Base::to_string() + ", " + std::to_string(m_offset) + ", " + std::to_string(m_size) +
             ")";
    }
 
  private:
    size_t m_offset;
    size_t m_size;
  };
 
  template <typename NT = double>
    requires std::is_floating_point_v<NT>
  class LinearCoordinates1D
  {
  public:
    using ctype = NT;
    static constexpr size_t dim = 1;
 
    LinearCoordinates1D(size_t grid_extent, double start, double stop)
        : start(start), stop(stop), grid_extent(grid_extent)
    {
      if (grid_extent == 0) throw std::runtime_error("LinearCoordinates1D: grid_extent must be > 0");
      a = (stop - start) / (grid_extent - 1.);
    }
 
    template <typename NT2>
    LinearCoordinates1D(const LinearCoordinates1D<NT2> &other)
        : LinearCoordinates1D(other.size(), other.start, other.stop)
    {
    }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION from_linear_index(size_t i) const
    {
      return device::array<size_t, 1>{i};
    }
 
    template <typename IT> NT KOKKOS_FORCEINLINE_FUNCTION forward(const IT &x) const { return start + a * x; }
 
    template <typename IT> device::array<NT, 1> KOKKOS_FORCEINLINE_FUNCTION forward(const device::array<IT, 1> &x) const
    {
      return {forward(x[0])};
    }
 
    NT KOKKOS_FORCEINLINE_FUNCTION backward(const NT &y) const { return (y - start) / a; }
 
    size_t KOKKOS_FORCEINLINE_FUNCTION size() const { return grid_extent; }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION sizes() const { return {grid_extent}; }
 
    const NT start, stop;
 
    template <typename NT2>
    friend bool operator==(const LinearCoordinates1D<NT> &lhs, const LinearCoordinates1D<NT2> &rhs)
    {
      if constexpr (!std::is_same_v<NT, NT2>) return false; // Different types, cannot be equal
      return lhs.start == rhs.start && lhs.stop == rhs.stop && lhs.grid_extent == rhs.grid_extent;
    }
 
    std::string to_string() const
    {
      return "LinearCoordinates1D(" + std::to_string(grid_extent) + ", " + std::to_string(start) + ", " +
             std::to_string(stop) + ")";
    }
 
  private:
    const size_t grid_extent;
    NT a;
  };
 
  template <typename NT = double>
    requires std::is_floating_point_v<NT>
  class LinearPeriodicCoordinates1D
  {
  public:
    using ctype = NT;
    static constexpr size_t dim = 1;
    static constexpr bool periodic = true;
 
    LinearPeriodicCoordinates1D(size_t grid_extent, double start, double stop)
        : start(start), stop(stop), grid_extent(grid_extent), extent(static_cast<NT>(grid_extent))
    {
      if (grid_extent == 0) throw std::runtime_error("LinearPeriodicCoordinates1D: grid_extent must be > 0");
      a = (stop - start) / NT(grid_extent);
    }
 
    template <typename NT2>
    LinearPeriodicCoordinates1D(const LinearPeriodicCoordinates1D<NT2> &other)
        : LinearPeriodicCoordinates1D(other.size(), other.start, other.stop)
    {
    }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION from_linear_index(size_t i) const
    {
      return device::array<size_t, 1>{i};
    }
 
    template <typename IT> NT KOKKOS_FORCEINLINE_FUNCTION forward(const IT &x) const { return start + a * x; }
 
    template <typename IT> device::array<NT, 1> KOKKOS_FORCEINLINE_FUNCTION forward(const device::array<IT, 1> &x) const
    {
      return {forward(x[0])};
    }
 
    NT KOKKOS_FORCEINLINE_FUNCTION backward(const NT &y) const
    {
      using Kokkos::floor;
      // fmod would keep the sign of idx, so fold explicitly into [0, grid_extent)
      NT idx = (y - start) / a;
      idx -= floor(idx / extent) * extent;
      // right at the seam the division can round such that the fold leaves idx just outside the interval
      if (idx >= extent) idx -= extent;
      if (idx < NT(0)) idx = NT(0);
      return idx;
    }
 
    size_t KOKKOS_FORCEINLINE_FUNCTION size() const { return grid_extent; }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION sizes() const { return {grid_extent}; }
 
    const NT start, stop;
 
    template <typename NT2>
    friend bool operator==(const LinearPeriodicCoordinates1D<NT> &lhs, const LinearPeriodicCoordinates1D<NT2> &rhs)
    {
      if constexpr (!std::is_same_v<NT, NT2>) return false; // Different types, cannot be equal
      return lhs.start == rhs.start && lhs.stop == rhs.stop && lhs.grid_extent == rhs.grid_extent;
    }
 
    std::string to_string() const
    {
      return "LinearPeriodicCoordinates1D(" + std::to_string(grid_extent) + ", " + std::to_string(start) + ", " +
             std::to_string(stop) + ")";
    }
 
  private:
    const size_t grid_extent;
    const NT extent;
    NT a;
  };
 
  template <typename NT = double>
    requires std::is_floating_point_v<NT>
  class LogarithmicCoordinates1D
  {
  public:
    using ctype = NT;
    static constexpr size_t dim = 1;
 
    LogarithmicCoordinates1D(size_t grid_extent, NT start, NT stop, NT bias)
        : start(start), stop(stop), bias(bias), grid_extent(grid_extent), gem1(grid_extent - 1.),
          gem1inv(1. / (grid_extent - 1.))
    {
      if (grid_extent == 0) throw std::runtime_error("LogarithmicCoordinates1D: grid_extent must be > 0");
      using Kokkos::expm1;
      a = bias;
      b = (stop - start) / expm1(a);
      c = start;
    }
 
    template <typename NT2>
    LogarithmicCoordinates1D(const LogarithmicCoordinates1D<NT2> &other)
        : LogarithmicCoordinates1D(other.size(), other.start, other.stop, other.bias)
    {
    }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION from_linear_index(size_t i) const
    {
      return device::array<size_t, 1>{i};
    }
 
    template <typename IT> NT KOKKOS_FORCEINLINE_FUNCTION forward(const IT &x) const
    {
      using Kokkos::expm1;
      return b * expm1(a * x * gem1inv) + c;
    }
 
    template <typename IT> device::array<NT, 1> KOKKOS_FORCEINLINE_FUNCTION forward(const device::array<IT, 1> &x) const
    {
      return {forward(x[0])};
    }
 
    NT KOKKOS_FORCEINLINE_FUNCTION backward(const NT &y) const
    {
      using Kokkos::log1p;
      return log1p((y - c) / b) * gem1 / a;
    }
 
    NT KOKKOS_FORCEINLINE_FUNCTION backward_derivative(const NT &y) const { return 1. / (y - c) * gem1 / a; }
 
    size_t KOKKOS_FORCEINLINE_FUNCTION size() const { return grid_extent; }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION sizes() const { return {grid_extent}; }
 
    const NT start, stop, bias;
 
    template <typename NT2>
    friend bool operator==(const LogarithmicCoordinates1D<NT> &lhs, const LogarithmicCoordinates1D<NT2> &rhs)
    {
      if constexpr (!std::is_same_v<NT, NT2>) return false; // Different types, cannot be equal
      return lhs.start == rhs.start && lhs.stop == rhs.stop && lhs.bias == rhs.bias &&
             lhs.grid_extent == rhs.grid_extent;
    }
 
    std::string to_string() const
    {
      return "LogarithmicCoordinates1D(" + std::to_string(start) + ", " + std::to_string(stop) + ", " +
             std::to_string(bias) + ", " + std::to_string(grid_extent) + ")";
    }
 
  private:
    const size_t grid_extent;
    const NT gem1, gem1inv;
    NT a, b, c;
  };
 
  template <typename NT = double>
    requires std::is_floating_point_v<NT>
  class FocusedLogCoordinates1D
  {
  public:
    using ctype = NT;
    static constexpr size_t dim = 1;
 
    FocusedLogCoordinates1D(size_t grid_extent, NT start, NT stop, NT center, NT focus)
        : start(start), stop(stop), center(center), focus(focus), grid_extent(grid_extent)
    {
      if (grid_extent < 2) throw std::runtime_error("FocusedLogCoordinates1D: grid_extent must be > 1");
      if (!(start > NT(0))) throw std::runtime_error("FocusedLogCoordinates1D: start must be > 0");
      if (!(stop > start)) throw std::runtime_error("FocusedLogCoordinates1D: stop must be > start");
      if (!(center > NT(0))) throw std::runtime_error("FocusedLogCoordinates1D: center must be > 0");
      if (!(focus >= NT(0))) throw std::runtime_error("FocusedLogCoordinates1D: focus must be >= 0");
 
      using Kokkos::log;
      const NT u_start = log(start);
      const NT u_stop = log(stop);
 
      u0 = log(center);
      c = focus;
      // sinh(s)/c is 0/0 at c == 0; below this threshold the map is a pure logarithmic grid to
      // within round-off anyway, so switch to it explicitly.
      pure_log = !(focus > NT(1e-8));
 
      if (pure_log) {
        // in this branch s *is* u, so forward() degenerates to start * exp(u - log(start))
        c_inv = NT(1);
        s_min = u_start;
        a = (u_stop - u_start) / NT(grid_extent - 1);
        g_min = s_min;
      } else {
        c_inv = NT(1) / c;
        s_min = s_of_u(u_start);
        a = (s_of_u(u_stop) - s_min) / NT(grid_extent - 1);
        using Kokkos::sinh;
        // Store sinh(s_min), NOT sinh(s_min) * c_inv: forward() subtracts this BEFORE scaling by
        // c_inv, which is what makes x == 0 cancel exactly. See the note in forward().
        g_min = sinh(s_min);
      }
      a_inv = NT(1) / a;
    }
 
    template <typename NT2>
    FocusedLogCoordinates1D(const FocusedLogCoordinates1D<NT2> &other)
        : FocusedLogCoordinates1D(other.size(), other.start, other.stop, other.center, other.focus)
    {
    }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION from_linear_index(size_t i) const
    {
      return device::array<size_t, 1>{i};
    }
 
    template <typename IT> NT KOKKOS_FORCEINLINE_FUNCTION forward(const IT &x) const
    {
      using Kokkos::exp, Kokkos::sinh;
      const NT s = s_min + a * x;
      // Factored as start * exp(...) rather than exp(log(start) + ...) so that x == 0 gives back
      // start bit-for-bit: the models pair grid element 0 with the p_grid_min parameter directly.
      // The upper end is then accurate to a handful of ulp, which nothing reads that way.
      //
      // The subtraction MUST happen before the multiplication by c_inv. Written the other way
      // round, `sinh(s) * c_inv - g_min` is a multiply-add, which the compiler is free to
      // contract into an FMA: the product is then kept at infinite precision and the ALREADY
      // ROUNDED g_min subtracted from it, so x == 0 yields that rounding error instead of zero
      // and start comes back off by an ulp. As written, `sinh(s) - g_min` is exactly zero at
      // x == 0 (identical argument, same function), and 0 * c_inv is exactly zero.
      const NT d = pure_log ? (s - g_min) : (sinh(s) - g_min) * c_inv;
      return start * exp(d);
    }
 
    template <typename IT> device::array<NT, 1> KOKKOS_FORCEINLINE_FUNCTION forward(const device::array<IT, 1> &x) const
    {
      return {forward(x[0])};
    }
 
    NT KOKKOS_FORCEINLINE_FUNCTION backward(const NT &y) const
    {
      using Kokkos::log;
      const NT u = log(y);
      return ((pure_log ? u : s_of_u(u)) - s_min) * a_inv;
    }
 
    NT KOKKOS_FORCEINLINE_FUNCTION backward_derivative(const NT &y) const
    {
      using Kokkos::log, Kokkos::sqrt;
      if (pure_log) return a_inv / y;
      const NT v = c * (log(y) - u0);
      return a_inv * c / (sqrt(v * v + NT(1)) * y);
    }
 
    size_t KOKKOS_FORCEINLINE_FUNCTION size() const { return grid_extent; }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION sizes() const { return {grid_extent}; }
 
    const NT start, stop, center, focus;
 
    template <typename NT2>
    friend bool operator==(const FocusedLogCoordinates1D<NT> &lhs, const FocusedLogCoordinates1D<NT2> &rhs)
    {
      if constexpr (!std::is_same_v<NT, NT2>) return false; // Different types, cannot be equal
      return lhs.start == rhs.start && lhs.stop == rhs.stop && lhs.center == rhs.center && lhs.focus == rhs.focus &&
             lhs.size() == rhs.size();
    }
 
    std::string to_string() const
    {
      return "FocusedLogCoordinates1D(" + std::to_string(grid_extent) + ", " + round_trip(start) + ", " +
             round_trip(stop) + ", " + round_trip(center) + ", " + round_trip(focus) + ")";
    }
 
  private:
    static std::string round_trip(const NT value)
    {
      char buffer[40];
      for (int precision = 6; precision < std::numeric_limits<NT>::max_digits10; ++precision) {
        std::snprintf(buffer, sizeof(buffer), "%.*g", precision, double(value));
        if (NT(std::strtod(buffer, nullptr)) == value) return std::string(buffer);
      }
      std::snprintf(buffer, sizeof(buffer), "%.*g", std::numeric_limits<NT>::max_digits10, double(value));
      return std::string(buffer);
    }
 
    NT KOKKOS_FORCEINLINE_FUNCTION s_of_u(const NT u) const
    {
      using Kokkos::log, Kokkos::sqrt;
      const NT v = c * (u - u0);
      const NT D = sqrt(v * v + NT(1));
      return v >= NT(0) ? log(v + D) : -log(D - v);
    }
 
    const size_t grid_extent;
    NT u0, s_min, a, a_inv, c, c_inv, g_min;
    bool pure_log;
  };
 
  template <typename Coordinates> auto make_grid(const Coordinates &coordinates)
  {
    using ctype = typename Coordinates::ctype;
    using cortype = device::array<ctype, Coordinates::dim>;
    std::vector<cortype> grid(coordinates.size());
    for (size_t i = 0; i < coordinates.size(); ++i) {
      const auto forwarded = coordinates.forward(coordinates.from_linear_index(i));
      for (size_t j = 0; j < Coordinates::dim; ++j) {
        grid[i][j] = forwarded[j];
      }
    }
    return grid;
  }
 
  template <typename Coordinates> auto make_idx_grid(const Coordinates &coordinates) -> std::vector<double>
  {
    using ctype = typename Coordinates::ctype;
    using cortype = device::array<ctype, Coordinates::dim>;
    std::vector<cortype> grid(coordinates.size());
    for (size_t i = 0; i < coordinates.size(); ++i) {
      const auto forwarded = coordinates.from_linear_index(i);
      for (size_t j = 0; j < Coordinates::dim; ++j) {
        grid[i][j] = forwarded[j];
      }
    }
    return grid;
  }
 
  template <typename Coordinates> std::vector<typename Coordinates::ctype> dump_grid(const Coordinates &coordinates)
  {
    using ctype = typename Coordinates::ctype;
    std::vector<typename Coordinates::ctype> grid(coordinates.size() * Coordinates::dim);
    for (size_t i = 0; i < coordinates.size(); ++i) {
      for (size_t j = 0; j < Coordinates::dim; ++j) {
        const auto forwarded = coordinates.forward(coordinates.from_linear_index(i));
        grid[i * coordinates.dim + j] = forwarded[j];
      }
    }
    return grid;
  }
 
  // Definitions of useful combined coordinates
  using LogCoordinates = LogarithmicCoordinates1D<double>;
  using LinCoordinates = LinearCoordinates1D<double>;
  using LogLogCoordinates = CoordinatePackND<LogarithmicCoordinates1D<double>, LogarithmicCoordinates1D<double>>;
  using LogLinCoordinates = CoordinatePackND<LogarithmicCoordinates1D<double>, LinearCoordinates1D<double>>;
  using LinLogCoordinates = CoordinatePackND<LinearCoordinates1D<double>, LogarithmicCoordinates1D<double>>;
  using LinLinCoordinates = CoordinatePackND<LinearCoordinates1D<double>, LinearCoordinates1D<double>>;
  using LogLogLinCoordinates =
      CoordinatePackND<LogarithmicCoordinates1D<double>, LogarithmicCoordinates1D<double>, LinearCoordinates1D<double>>;
  using LogLinLinCoordinates =
      CoordinatePackND<LogarithmicCoordinates1D<double>, LinearCoordinates1D<double>, LinearCoordinates1D<double>>;
  using LinLinLinCoordinates =
      CoordinatePackND<LinearCoordinates1D<double>, LinearCoordinates1D<double>, LinearCoordinates1D<double>>;
 
  // Logarithmic coordinates focused on an interior scale, see FocusedLogCoordinates1D
  using FocusedLogCoordinates = FocusedLogCoordinates1D<double>;
  using FocusedLogLinCoordinates = CoordinatePackND<FocusedLogCoordinates1D<double>, LinearCoordinates1D<double>>;
  using FocusedLogLinLinCoordinates = CoordinatePackND<FocusedLogCoordinates1D<double>, LinearCoordinates1D<double>,
                                                       LinearCoordinates1D<double>>;
 
  // Periodic (angular) coordinates
  using LinPeriodicCoordinates = LinearPeriodicCoordinates1D<double>;
  using LogLinPeriodicCoordinates =
      CoordinatePackND<LogarithmicCoordinates1D<double>, LinearPeriodicCoordinates1D<double>>;
  using LogLinLinPeriodicCoordinates = CoordinatePackND<LogarithmicCoordinates1D<double>, LinearCoordinates1D<double>,
                                                        LinearPeriodicCoordinates1D<double>>;
  using FocusedLogLinPeriodicCoordinates =
      CoordinatePackND<FocusedLogCoordinates1D<double>, LinearPeriodicCoordinates1D<double>>;
  using FocusedLogLinLinPeriodicCoordinates =
      CoordinatePackND<FocusedLogCoordinates1D<double>, LinearCoordinates1D<double>,
                       LinearPeriodicCoordinates1D<double>>;
} // namespace DiFfRG