diff --git a/example/poisson/poisson_package.cpp b/example/poisson/poisson_package.cpp
index e3f5bf75d3ae..00cd5071af37 100644
--- a/example/poisson/poisson_package.cpp
+++ b/example/poisson/poisson_package.cpp
@@ -160,9 +160,9 @@ TaskStatus SumMass(T *u, Real *reduce_sum) {
 
   const parthenon::Coordinates_t &coords = GetCoords(pm);
   const int ndim = v.GetNdim();
-  const Real dx = coords.Dxc<X1DIR>();
+  const Real dx = coords.Dxc<X1DIR>(0, 0, 0);
   for (int i = X2DIR; i <= ndim; i++) {
-    const Real dy = coords.DxcFA(i);
+    const Real dy = coords.Dxc(i, 0, 0, 0);
     PARTHENON_REQUIRE_THROWS(dx == dy,
                              "SumMass requires that DX be equal in all directions.");
   }
@@ -231,9 +231,9 @@ TaskStatus UpdatePhi(T *u, T *du) {
 
   const parthenon::Coordinates_t &coords = GetCoords(pm);
   const int ndim = v.GetNdim();
-  const Real dx = coords.Dxc<X1DIR>();
+  const Real dx = coords.Dxc<X1DIR>(0, 0, 0);
   for (int i = X2DIR; i <= ndim; i++) {
-    const Real dy = coords.DxcFA(i);
+    const Real dy = coords.Dxc(i, 0, 0, 0);
     PARTHENON_REQUIRE_THROWS(dx == dy,
                              "UpdatePhi requires that DX be equal in all directions.");
   }
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 1ea2e8bcd0bc..7295c3adc388 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -90,7 +90,10 @@ set(COMPILED_WITH ${CMAKE_CXX_COMPILER})
 set(COMPILER_COMMAND "<not-implemented>") # TODO: Put something more descriptive here
 set(COMPILER_FLAGS "<not-implemented>") # TODO: Put something more descriptive here
 
-set(COORDINATE_TYPE UniformCartesian) # TODO: Make this an option when more are available
+if (NOT COORDINATE_TYPE)
+  set(COORDINATE_TYPE UniformCartesian)
+endif()
+message(STATUS "COORDINATE_TYPE = " ${COORDINATE_TYPE})
 
 configure_file(config.hpp.in generated/config.hpp @ONLY)
 
@@ -115,6 +118,8 @@ add_library(parthenon
 
   coordinates/coordinates.hpp
   coordinates/uniform_cartesian.hpp
+  coordinates/uniform_coordinates.hpp
+  coordinates/uniform_spherical.hpp
 
   driver/driver.cpp
   driver/driver.hpp
diff --git a/src/coordinates/coordinates.hpp b/src/coordinates/coordinates.hpp
index d1290deea425..e89a008f8da0 100644
--- a/src/coordinates/coordinates.hpp
+++ b/src/coordinates/coordinates.hpp
@@ -16,6 +16,7 @@
 #include "config.hpp"
 
 #include "uniform_cartesian.hpp"
+#include "uniform_spherical.hpp"
 
 namespace parthenon {
 
diff --git a/src/coordinates/uniform_cartesian.hpp b/src/coordinates/uniform_cartesian.hpp
index 04cec9d41974..32f38bd6a9b5 100644
--- a/src/coordinates/uniform_cartesian.hpp
+++ b/src/coordinates/uniform_cartesian.hpp
@@ -13,319 +13,20 @@
 #ifndef COORDINATES_UNIFORM_CARTESIAN_HPP_
 #define COORDINATES_UNIFORM_CARTESIAN_HPP_
 
-#include <array>
-#include <cassert>
-
-#include "basic_types.hpp"
-#include "defs.hpp"
-#include "globals.hpp"
-
-#include <Kokkos_Macros.hpp>
+#include "uniform_coordinates.hpp"
 
 namespace parthenon {
 
-class UniformCartesian {
+class UniformCartesian : public UniformCoordinates<UniformCartesian> {
  public:
   UniformCartesian() = default;
-  UniformCartesian(const RegionSize &rs, ParameterInput *pin) {
-    for (auto &dir : {X1DIR, X2DIR, X3DIR}) {
-      dx_[dir - 1] = (rs.xmax(dir) - rs.xmin(dir)) / rs.nx(dir);
-      istart_[dir - 1] = (!rs.symmetry(dir) ? Globals::nghost : 0);
-      xmin_[dir - 1] = rs.xmin(dir) - istart_[dir - 1] * dx_[dir - 1];
-    }
-    area_[0] = dx_[1] * dx_[2];
-    area_[1] = dx_[0] * dx_[2];
-    area_[2] = dx_[0] * dx_[1];
-    cell_volume_ = dx_[0] * dx_[1] * dx_[2];
-  }
+  UniformCartesian(const RegionSize &rs, ParameterInput *pin) 
+    : UniformCoordinates<UniformCartesian>(rs, pin) {}
   UniformCartesian(const UniformCartesian &src, int coarsen)
-      : istart_(src.GetStartIndex()) {
-    dx_ = src.Dx_();
-    xmin_ = src.GetXmin();
-    xmin_[0] += istart_[0] * dx_[0] * (1 - coarsen);
-    xmin_[1] += istart_[1] * dx_[1] * (1 - coarsen);
-    xmin_[2] += istart_[2] * dx_[2] * (1 - coarsen);
-    dx_[0] *= coarsen;
-    dx_[1] *= (istart_[1] > 0 ? coarsen : 1);
-    dx_[2] *= (istart_[2] > 0 ? coarsen : 1);
-    area_[0] = dx_[1] * dx_[2];
-    area_[1] = dx_[0] * dx_[2];
-    area_[2] = dx_[0] * dx_[1];
-    cell_volume_ = dx_[0] * dx_[1] * dx_[2];
-  }
-
-  //----------------------------------------
-  // Dxc: Distance between cell centers
-  //----------------------------------------
-  template <int dir>
-  KOKKOS_FORCEINLINE_FUNCTION Real Dxc() const {
-    assert(dir > 0 && dir < 4);
-    return dx_[dir - 1];
-  }
-  template <int dir, class... Args>
-  KOKKOS_FORCEINLINE_FUNCTION Real Dxc(Args... args) const {
-    assert(dir > 0 && dir < 4);
-    return dx_[dir - 1];
-  }
-  template <class... Args>
-  KOKKOS_FORCEINLINE_FUNCTION Real DxcFA(const int dir, Args... args) const {
-    assert(dir > 0 && dir < 4);
-    return dx_[dir - 1];
-  }
-
-  //----------------------------------------
-  // Dxf: Distance between cell faces
-  //----------------------------------------
-  template <int dir, int face, class... Args>
-  KOKKOS_FORCEINLINE_FUNCTION Real Dxf(Args... args) const {
-    assert(dir > 0 && dir < 4 && face > 0 && face < 4);
-    return dx_[face - 1];
-  }
-  template <int dir, class... Args>
-  KOKKOS_FORCEINLINE_FUNCTION Real Dxf(Args... args) const {
-    assert(dir > 0 && dir < 4);
-    return dx_[dir - 1];
-  }
-
-  //----------------------------------------
-  // Xf: Positions at cell centers
-  //----------------------------------------
-  template <int dir>
-  KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int idx) const {
-    assert(dir > 0 && dir < 4);
-    return xmin_[dir - 1] + (idx + 0.5) * dx_[dir - 1];
-  }
-  template <int dir>
-  KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int k, const int j, const int i) const {
-    assert(dir > 0 && dir < 4);
-    switch (dir) {
-    case 1:
-      return Xc<dir>(i);
-    case 2:
-      return Xc<dir>(j);
-    case 3:
-      return Xc<dir>(k);
-    default:
-      PARTHENON_FAIL("Unknown dir");
-      return 0; // To appease compiler
-    }
-  }
-
-  //----------------------------------------
-  // Xf: Positions on Faces
-  //----------------------------------------
-  template <int dir, int face>
-  KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int idx) const {
-    assert(dir > 0 && dir < 4 && face > 0 && face < 4);
-    // Return position in direction "dir" along index "idx" on face "face"
-    if constexpr (dir == face) {
-      return xmin_[dir - 1] + idx * dx_[dir - 1];
-    } else {
-      return Xc<dir>(idx);
-    }
-  }
-  template <int dir, int face>
-  KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int k, const int j, const int i) const {
-    assert(dir > 0 && dir < 4);
-    switch (dir) {
-    case 1:
-      return Xf<dir, face>(i);
-    case 2:
-      return Xf<dir, face>(j);
-    case 3:
-      return Xf<dir, face>(k);
-    default:
-      PARTHENON_FAIL("Unknown dir");
-      return 0; // To appease compiler
-    }
-  }
-
-  template <int dir>
-  KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int idx) const {
-    assert(dir > 0 && dir < 4);
-    // Return position in direction "dir" along index "idx" on face "dir"
-    return xmin_[dir - 1] + idx * dx_[dir - 1];
-  }
-  template <int dir>
-  KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int k, const int j, const int i) const {
-    assert(dir > 0 && dir < 4);
-    switch (dir) {
-    case 1:
-      return Xf<dir>(i);
-    case 2:
-      return Xf<dir>(j);
-    case 3:
-      return Xf<dir>(k);
-    default:
-      PARTHENON_FAIL("Unknown dir");
-      return 0; // To appease compiler
-    }
-  }
-
-  template <int dir, TopologicalElement el>
-  KOKKOS_FORCEINLINE_FUNCTION Real X(const int idx) const {
-    if constexpr (dir == X1DIR && TopologicalOffsetI(el)) {
-      return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2
-    } else if constexpr (dir == X2DIR && TopologicalOffsetJ(el)) {
-      return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2
-    } else if constexpr (dir == X3DIR && TopologicalOffsetK(el)) {
-      return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2
-    } else {
-      return xmin_[dir - 1] + (idx + 0.5) * dx_[dir - 1]; // idx
-    }
-    return 0; // This should never be reached, but w/o it some compilers generate warnings
-  }
-
-  template <int dir, TopologicalElement el>
-  KOKKOS_FORCEINLINE_FUNCTION Real X(const int k, const int j, const int i) const {
-    assert(dir > 0 && dir < 4);
-    switch (dir) {
-    case 1:
-      return X<dir, el>(i);
-    case 2:
-      return X<dir, el>(j);
-    case 3:
-      return X<dir, el>(k);
-    default:
-      PARTHENON_FAIL("Unknown dir");
-      return 0; // To appease compiler
-    }
-  }
-
-  //----------------------------------------
-  // CellWidth: Width of cells at cell centers
-  //----------------------------------------
-  template <int dir, class... Args>
-  KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(Args... args) const {
-    assert(dir > 0 && dir < 4);
-    return dx_[dir - 1];
-  }
-  template <class... Args>
-  KOKKOS_FORCEINLINE_FUNCTION Real CellWidthFA(const int dir, Args... args) const {
-    assert(dir > 0 && dir < 4);
-    return dx_[dir - 1];
-  }
+      : UniformCoordinates<UniformCartesian>(src, coarsen) {}
 
-  //----------------------------------------
-  // EdgeLength: Length of cell edges
-  //----------------------------------------
-  template <int dir, class... Args>
-  KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(Args... args) const {
-    assert(dir > 0 && dir < 4);
-    return CellWidth<dir>();
-  }
-  template <class... Args>
-  KOKKOS_FORCEINLINE_FUNCTION Real EdgeLengthFA(const int dir, Args... args) const {
-    return CellWidthFA(dir);
-  }
-
-  //----------------------------------------
-  // FaceArea: Area of cell areas
-  //----------------------------------------
-  template <int dir, class... Args>
-  KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(Args... args) const {
-    assert(dir > 0 && dir < 4);
-    return area_[dir - 1];
-  }
-  template <class... Args>
-  KOKKOS_FORCEINLINE_FUNCTION Real FaceAreaFA(const int dir, Args... args) const {
-    assert(dir > 0 && dir < 4);
-    return area_[dir - 1];
-  }
-
-  //----------------------------------------
-  // CellVolume
-  //----------------------------------------
-  template <class... Args>
-  KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(Args... args) const {
-    return cell_volume_;
-  }
-
-  //----------------------------------------
-  // Generalized volume
-  //----------------------------------------
-  template <TopologicalElement el, class... Args>
-  KOKKOS_FORCEINLINE_FUNCTION Real Volume(Args... args) const {
-    using TE = TopologicalElement;
-    if constexpr (el == TE::CC) {
-      return cell_volume_;
-    } else if constexpr (el == TE::F1) {
-      return area_[X1DIR - 1];
-    } else if constexpr (el == TE::F2) {
-      return area_[X2DIR - 1];
-    } else if constexpr (el == TE::F3) {
-      return area_[X3DIR - 1];
-    } else if constexpr (el == TE::E1) {
-      return dx_[X1DIR - 1];
-    } else if constexpr (el == TE::E2) {
-      return dx_[X2DIR - 1];
-    } else if constexpr (el == TE::E3) {
-      return dx_[X3DIR - 1];
-    } else if constexpr (el == TE::NN) {
-      return 1.0;
-    }
-    PARTHENON_FAIL("If you reach this point, someone has added a new value to the the "
-                   "TopologicalElement enum.");
-    return 0.0;
-  }
-
-  template <class... Args>
-  KOKKOS_FORCEINLINE_FUNCTION Real Volume(CellLevel cl, TopologicalElement el,
-                                          Args... args) const {
-    using TE = TopologicalElement;
-    if (cl == CellLevel::same) {
-      if (el == TE::CC) {
-        return cell_volume_;
-      } else if (el == TE::F1) {
-        return area_[X1DIR - 1];
-      } else if (el == TE::F2) {
-        return area_[X2DIR - 1];
-      } else if (el == TE::F3) {
-        return area_[X3DIR - 1];
-      } else if (el == TE::E1) {
-        return dx_[X1DIR - 1];
-      } else if (el == TE::E2) {
-        return dx_[X2DIR - 1];
-      } else if (el == TE::E3) {
-        return dx_[X3DIR - 1];
-      } else if (el == TE::NN) {
-        return 1.0;
-      }
-    } else if (cl == CellLevel::fine) {
-      if (el == TE::CC) {
-        return cell_volume_ / 8.0;
-      } else if (el == TE::F1) {
-        return area_[X1DIR - 1] / 4.0;
-      } else if (el == TE::F2) {
-        return area_[X2DIR - 1] / 4.0;
-      } else if (el == TE::F3) {
-        return area_[X3DIR - 1] / 4.0;
-      } else if (el == TE::E1) {
-        return dx_[X1DIR - 1] / 2.0;
-      } else if (el == TE::E2) {
-        return dx_[X2DIR - 1] / 2.0;
-      } else if (el == TE::E3) {
-        return dx_[X3DIR - 1] / 2.0;
-      } else if (el == TE::NN) {
-        return 1.0;
-      }
-    }
-    PARTHENON_FAIL("If you reach this point, someone has added a new value to the the "
-                   "TopologicalElement enum.");
-    return 0.0;
-  }
-
-  const std::array<Real, 3> &GetXmin() const { return xmin_; }
-  const std::array<int, 3> &GetStartIndex() const { return istart_; }
-  const char *Name() const { return name_; }
-
- private:
-  std::array<int, 3> istart_;
-  std::array<Real, 3> xmin_, dx_, area_;
-  Real cell_volume_;
   constexpr static const char *name_ = "UniformCartesian";
-
-  const std::array<Real, 3> &Dx_() const { return dx_; }
+ private:
 };
 
 } // namespace parthenon
diff --git a/src/coordinates/uniform_coordinates.hpp b/src/coordinates/uniform_coordinates.hpp
new file mode 100644
index 000000000000..93dd6a1f5ebe
--- /dev/null
+++ b/src/coordinates/uniform_coordinates.hpp
@@ -0,0 +1,343 @@
+//========================================================================================
+// (C) (or copyright) 2024. 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 COORDINATES_UNIFORM_COORDINATES_HPP_
+#define COORDINATES_UNIFORM_COORDINATES_HPP_
+
+#include <array>
+#include <cassert>
+
+#include "basic_types.hpp"
+#include "defs.hpp"
+#include "globals.hpp"
+
+#include <Kokkos_Macros.hpp>
+
+namespace parthenon {
+
+template <typename System>
+class UniformCoordinates {
+ public:
+  UniformCoordinates() = default;
+  UniformCoordinates(const RegionSize &rs, ParameterInput *pin) {
+    for (auto &dir : {X1DIR, X2DIR, X3DIR}) {
+      dx_[dir - 1] = (rs.xmax(dir) - rs.xmin(dir)) / rs.nx(dir);
+      istart_[dir - 1] = (!rs.symmetry(dir) ? Globals::nghost : 0);
+      xmin_[dir - 1] = rs.xmin(dir) - istart_[dir - 1] * dx_[dir - 1];
+    }
+    area_[0] = dx_[1] * dx_[2];
+    area_[1] = dx_[0] * dx_[2];
+    area_[2] = dx_[0] * dx_[1];
+    cell_volume_ = dx_[0] * dx_[1] * dx_[2];
+  }
+  UniformCoordinates(const UniformCoordinates<System> &src, int coarsen)
+      : istart_(src.GetStartIndex()) {
+    dx_ = src.Dx_();
+    xmin_ = src.GetXmin();
+    xmin_[0] += istart_[0] * dx_[0] * (1 - coarsen);
+    xmin_[1] += istart_[1] * dx_[1] * (1 - coarsen);
+    xmin_[2] += istart_[2] * dx_[2] * (1 - coarsen);
+    dx_[0] *= coarsen;
+    dx_[1] *= (istart_[1] > 0 ? coarsen : 1);
+    dx_[2] *= (istart_[2] > 0 ? coarsen : 1);
+    area_[0] = dx_[1] * dx_[2];
+    area_[1] = dx_[0] * dx_[2];
+    area_[2] = dx_[0] * dx_[1];
+    cell_volume_ = dx_[0] * dx_[1] * dx_[2];
+  }
+  //----------------------------------------
+  // Dxc: Distance between cell centers
+  //----------------------------------------
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real Dx() const {
+    static_assert(dir > 0 && dir < 4);
+    return dx_[dir - 1];
+  }
+  KOKKOS_FORCEINLINE_FUNCTION Real Dx(const int dir) const {
+    assert(dir > 0 && dir < 4);
+    return dx_[dir - 1];
+  }
+
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int idx) const {
+    static_assert(dir > 0 && dir < 4);
+    return dx_[dir - 1];
+  }
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int k, const int j, const int i) const {
+    static_assert(dir > 0 && dir < 4);
+    if constexpr (dir == X1DIR) return static_cast<const System *>(this)->template Dxc<dir>(i);
+    else if constexpr (dir == X2DIR) return static_cast<const System *>(this)->template Dxc<dir>(j);
+    else if constexpr (dir == X3DIR) return static_cast<const System *>(this)->template Dxc<dir>(k);
+    PARTHENON_FAIL("Unknown dir.");
+    return 0.0;
+  }
+  KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int dir, const int k, const int j, const int i) const {
+    assert(dir > 0 && dir < 4);
+    if (dir == X1DIR) return static_cast<const System *>(this)->template Dxc<X1DIR>(k, j, i);
+    else if (dir == X2DIR) return static_cast<const System *>(this)->template Dxc<X2DIR>(k, j, i);
+    else if (dir == X3DIR) return static_cast<const System *>(this)->template Dxc<X3DIR>(k, j, i);
+    PARTHENON_FAIL("Unknown dir.");
+    return 0.0;
+  }
+
+  //----------------------------------------
+  // Dxf: Distance between cell faces
+  //----------------------------------------
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real Dxf(const int idx) const {
+    static_assert(dir > 0 && dir < 4);
+    return dx_[dir - 1];
+  }
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real Dxf(const int k, const int j, const int i) const {
+    static_assert(dir > 0 && dir < 4);
+    return dx_[dir - 1];
+  }
+
+  //----------------------------------------
+  // Xc: Positions at cell centroids
+  //----------------------------------------
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int idx) const {
+    static_assert(dir > 0 && dir < 4);
+    return xmin_[dir - 1] + (idx + 0.5) * dx_[dir - 1];
+  }
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int k, const int j, const int i) const {
+    static_assert(dir > 0 && dir < 4);
+    if constexpr (dir == X1DIR) return static_cast<const System *>(this)->template Xc<dir>(i);
+    else if constexpr (dir == X2DIR) return static_cast<const System *>(this)->template Xc<dir>(j);
+    else if constexpr (dir == X3DIR) return static_cast<const System *>(this)->template Xc<dir>(k);
+    return 0; // To appease compiler
+  }
+
+  //----------------------------------------
+  // Xf: Positions on Faces
+  //----------------------------------------
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int idx) const {
+    static_assert(dir > 0 && dir < 4);
+    return xmin_[dir - 1] + idx * dx_[dir - 1];
+  }
+  template <int dir, int face>
+  KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int k, const int j, const int i) const {
+    static_assert(dir > 0 && dir < 4);
+    static_assert(face > 0 && face < 4);
+    if constexpr (dir == X1DIR && dir == face) {
+      return Xf<dir>(i);
+    } else if constexpr (dir == X1DIR) {
+      return static_cast<const System *>(this)->template Xc<dir>(k, j, i);
+    } else if constexpr (dir == X2DIR && dir == face) {
+      return Xf<dir>(j);
+    } else if constexpr (dir == X2DIR) {
+      return static_cast<const System *>(this)->template Xc<dir>(k, j, i);
+    } else if constexpr (dir == X3DIR && dir == face) {
+      return Xf<dir>(k);
+    } else if constexpr (dir == X3DIR) {
+      return static_cast<const System *>(this)->template Xc<dir>(k, j, i);
+    }
+    return 0; // To appease compiler
+  }
+
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int k, const int j, const int i) const {
+    return Xf<dir, dir>(k, j, i);
+  }
+
+  template <int dir, TopologicalElement el>
+  KOKKOS_FORCEINLINE_FUNCTION Real X(const int idx) const {
+    static_assert(dir > 0 && dir < 4);
+    if constexpr (dir == X1DIR && TopologicalOffsetI(el)) {
+      return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2
+    } else if constexpr (dir == X2DIR && TopologicalOffsetJ(el)) {
+      return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2
+    } else if constexpr (dir == X3DIR && TopologicalOffsetK(el)) {
+      return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2
+    } else {
+      return static_cast<const System *>(this)->template Xc<dir>(idx); // idx
+    }
+    return 0; // This should never be reached, but w/o it some compilers generate warnings
+  }
+  template <int dir, TopologicalElement el>
+  KOKKOS_FORCEINLINE_FUNCTION Real X(const int k, const int j, const int i) const {
+    static_assert(dir > 0 && dir < 4);
+    if constexpr (dir == X1DIR) return X<dir, el>(i);
+    else if constexpr (dir == X2DIR) return X<dir, el>(j);
+    else if constexpr (dir == X3DIR) return X<dir, el>(k);
+    return 0.0;
+  }
+
+  template <int dir, TopologicalElement el>
+  KOKKOS_FORCEINLINE_FUNCTION
+  Real Scale(const int k, const int j, const int i) const {
+    if constexpr (dir > 0 && dir < 4) return 1.0;
+    PARTHENON_FAIL("Unknown dir");
+    return 0.0;
+  }
+
+  template <TopologicalElement el>
+  KOKKOS_FORCEINLINE_FUNCTION
+  Real Scale(const int dir, const int k, const int j, const int i) const {
+    assert(dir > 0 && dir < 4);
+    if (dir == X1DIR) return static_cast<const System *>(this)->template Scale<X1DIR, el>(k, j, i);
+    else if (dir == X2DIR) return static_cast<const System *>(this)->template Scale<X2DIR, el>(k, j, i);
+    else if (dir == X3DIR) return static_cast<const System *>(this)->template Scale<X3DIR, el>(k, j, i);
+    return 0.0;
+  }
+
+  //----------------------------------------
+  // CellWidth: Width of cells at cell centers
+  //----------------------------------------
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(const int k, const int j, const int i) const {
+    assert(dir > 0 && dir < 4);
+    return dx_[dir - 1];
+  }
+  KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(const int dir, const int k, const int j, const int i) const {
+    assert(dir > 0 && dir < 4);
+    if (dir == X1DIR) return static_cast<const System *>(this)->template CellWidth<X1DIR>(k, j, i);
+    else if (dir == X2DIR) return static_cast<const System *>(this)->template CellWidth<X2DIR>(k, j, i);
+    else if (dir == X3DIR) return static_cast<const System *>(this)->template CellWidth<X3DIR>(k, j, i);
+    return 0.0;
+  }
+
+  //----------------------------------------
+  // EdgeLength: Length of cell edges
+  //----------------------------------------
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int k, const int j, const int i) const {
+    assert(dir > 0 && dir < 4);
+    return CellWidth<dir>(k, j, i);
+  }
+  KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int dir, const int k, const int j, const int i) const {
+    return CellWidth(dir, k, j, i);
+  }
+
+  //----------------------------------------
+  // FaceArea: Area of cell areas
+  //----------------------------------------
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const int k, const int j, const int i) const {
+    static_assert(dir > 0 && dir < 4);
+    return area_[dir - 1];
+  }
+  KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const int dir, const int k, const int j, const int i) const {
+    assert(dir > 0 && dir < 4);
+    switch(dir) {
+    case X1DIR:
+      return static_cast<const System *>(this)->template FaceArea<X1DIR>(k, j, i);
+    case X2DIR:
+      return static_cast<const System *>(this)->template FaceArea<X2DIR>(k, j, i);
+    case X3DIR:
+      return static_cast<const System *>(this)->template FaceArea<X3DIR>(k, j, i);
+    default:
+      PARTHENON_FAIL("Unknown dir");
+      return 0; // To appease compiler
+    }
+  }
+
+  //----------------------------------------
+  // CellVolume
+  //----------------------------------------
+  KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, const int i) const {
+    return cell_volume_;
+  }
+
+  //----------------------------------------
+  // Generalized volume
+  //----------------------------------------
+  template <TopologicalElement el>
+  KOKKOS_FORCEINLINE_FUNCTION Real Volume(const int k, const int j, const int i) const {
+    using TE = TopologicalElement;
+    if constexpr (el == TE::CC) {
+      return static_cast<const System *>(this)->CellVolume(k, j, i);
+    } else if constexpr (el == TE::F1) {
+      return static_cast<const System *>(this)->template FaceArea<X1DIR>(k, j, i);
+    } else if constexpr (el == TE::F2) {
+      return static_cast<const System *>(this)->template FaceArea<X2DIR>(k, j, i);
+    } else if constexpr (el == TE::F3) {
+      return static_cast<const System *>(this)->template FaceArea<X3DIR>(k, j, i);
+    } else if constexpr (el == TE::E1) {
+      return static_cast<const System *>(this)->template EdgeLength<X1DIR>(k, j, i);
+    } else if constexpr (el == TE::E2) {
+      return static_cast<const System *>(this)->template EdgeLength<X2DIR>(k, j, i);
+    } else if constexpr (el == TE::E3) {
+      return static_cast<const System *>(this)->template EdgeLength<X3DIR>(k, j, i);
+    } else if constexpr (el == TE::NN) {
+      return 1.0;
+    }
+    PARTHENON_FAIL("If you reach this point, someone has added a new value to the the "
+                   "TopologicalElement enum.");
+    return 0.0;
+  }
+
+  template <class... Args>
+  KOKKOS_FORCEINLINE_FUNCTION Real Volume(CellLevel cl, TopologicalElement el,
+                                          const int k, const int j, const int i) const {
+    using TE = TopologicalElement;
+    if (cl == CellLevel::same) {
+      if (el == TE::CC) {
+        return cell_volume_;
+      } else if (el == TE::F1) {
+        return area_[X1DIR - 1];
+      } else if (el == TE::F2) {
+        return area_[X2DIR - 1];
+      } else if (el == TE::F3) {
+        return area_[X3DIR - 1];
+      } else if (el == TE::E1) {
+        return dx_[X1DIR - 1];
+      } else if (el == TE::E2) {
+        return dx_[X2DIR - 1];
+      } else if (el == TE::E3) {
+        return dx_[X3DIR - 1];
+      } else if (el == TE::NN) {
+        return 1.0;
+      }
+    } else if (cl == CellLevel::fine) {
+      if (el == TE::CC) {
+        return cell_volume_ / 8.0;
+      } else if (el == TE::F1) {
+        return area_[X1DIR - 1] / 4.0;
+      } else if (el == TE::F2) {
+        return area_[X2DIR - 1] / 4.0;
+      } else if (el == TE::F3) {
+        return area_[X3DIR - 1] / 4.0;
+      } else if (el == TE::E1) {
+        return dx_[X1DIR - 1] / 2.0;
+      } else if (el == TE::E2) {
+        return dx_[X2DIR - 1] / 2.0;
+      } else if (el == TE::E3) {
+        return dx_[X3DIR - 1] / 2.0;
+      } else if (el == TE::NN) {
+        return 1.0;
+      }
+    }
+    PARTHENON_FAIL("If you reach this point, someone has added a new value to the the "
+                   "TopologicalElement enum.");
+    return 0.0;
+  }
+
+  const std::array<Real, 3> &GetXmin() const { return xmin_; }
+  const std::array<int, 3> &GetStartIndex() const { return istart_; }
+  const char *Name() const { return System::name_; }
+ private:
+  std::array<int, 3> istart_;
+  std::array<Real, 3> xmin_, dx_, area_;
+  Real cell_volume_;
+
+  const std::array<Real, 3> &Dx_() const { return dx_; }
+};
+
+} // namespace parthenon
+
+#endif // COORDINATES_UNIFORM_COORDINATES_HPP_
diff --git a/src/coordinates/uniform_spherical.hpp b/src/coordinates/uniform_spherical.hpp
new file mode 100644
index 000000000000..6f8cae6dc938
--- /dev/null
+++ b/src/coordinates/uniform_spherical.hpp
@@ -0,0 +1,174 @@
+//========================================================================================
+// (C) (or copyright) 2023-2024. 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 COORDINATES_UNIFORM_SPHERICAL_HPP_
+#define COORDINATES_UNIFORM_SPHERICAL_HPP_
+
+#include "uniform_coordinates.hpp"
+
+namespace parthenon {
+
+class UniformSpherical : public UniformCoordinates<UniformSpherical> {
+ using base_t = UniformCoordinates<UniformSpherical>;
+ public:
+  using base_t::Dxc;
+  using base_t::Xc;
+  using base_t::Scale;
+  using base_t::CellWidth;
+  using base_t::Volume;
+  UniformSpherical() = default;
+  UniformSpherical(const RegionSize &rs, ParameterInput *pin) 
+    : UniformCoordinates<UniformSpherical>(rs, pin) {}
+  UniformSpherical(const UniformSpherical &src, int coarsen)
+      : UniformCoordinates<UniformSpherical>(src, coarsen) {}
+  constexpr static const char *name_ = "UniformSpherical";
+
+  //----------------------------------------
+  // Dxc: Distance between cell centers
+  //----------------------------------------
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int idx) const {
+    static_assert(dir > 0 && dir < 4);
+    return Xc<dir>(idx) - Xc<dir>(idx-1);
+  }
+
+  //----------------------------------------
+  // Xc: Positions at cell centroids
+  //----------------------------------------
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int idx) const {
+    static_assert(dir > 0 && dir < 4);
+    if constexpr (dir == X1DIR) {
+      const Real r0 = Xf<X1DIR>(idx);
+      const Real r0sq = r0 * r0;
+      const Real r1 = Xf<X1DIR>(idx + 1);
+      const Real r1sq = r1 * r1;
+      return 0.75 * (r1sq * r1sq - r0sq * r0sq) / (r1sq * r1 - r0sq * r0);
+    } else if constexpr (dir == X2DIR) {
+      const Real th0 = Xf<X2DIR>(idx);
+      const Real sth0 = std::sin(th0);
+      const Real cth0 = std::cos(th0);
+      const Real th1 = Xf<X2DIR>(idx + 1);
+      const Real sth1 = std::sin(th1);
+      const Real cth1 = std::cos(th1);
+      return (th0 * cth0 - th1 * cth1 - sth0 + sth1) / (cth0 - cth1);
+    }
+    return base_t::Xc<X3DIR>(idx);
+  }
+
+  template <int dir, TopologicalElement el>
+  KOKKOS_FORCEINLINE_FUNCTION
+  Real Scale(const int k, const int j, const int i) const {
+    static_assert(dir > 0 && dir < 4);
+    using TE = TopologicalElement;
+    if constexpr (dir == X1DIR) return 1.0;
+    else {
+      const Real r = X<X1DIR, el>(k, j, i);
+      if constexpr (dir == X2DIR) return r;
+      else {
+        const Real th = X<X2DIR, el>(k, j, i);
+        return r * std::sin(th);
+      }
+    }
+    return 0.0;
+  }
+
+  //----------------------------------------
+  // CellWidth: width of cell through the centroid
+  //----------------------------------------
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(const int k, const int j, const int i) const {
+    using TE = TopologicalElement;
+    static_assert(dir > 0 && dir < 4);
+    if constexpr (dir == X1DIR) return Dx<dir>();
+    else if constexpr (dir == X2DIR) return Xc<X1DIR>(i) * Dx<dir>();
+    else if constexpr (dir == X3DIR) return Xc<X1DIR>(i) * std::sin(Xc<X2DIR>(j)) * Dx<dir>();
+    return 0.0;
+  }
+
+  //----------------------------------------
+  // EdgeLength: Length of cell edges
+  //----------------------------------------
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int k, const int j, const int i) const {
+    static_assert(dir > 0 && dir < 4);
+    if constexpr (dir == X1DIR) {
+      // radial direction is trivial
+      return Dx<dir>();
+    } else if constexpr (dir == X2DIR) {
+      // theta direction
+      return Xf<X1DIR>(k, j, i) * Dx<dir>();
+    }
+    // phi direction
+    return Xf<X1DIR>(k, j, i) * std::sin(Xf<X2DIR>(k, j, i)) * Dx<dir>();
+  }
+  KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int dir, const int k, const int j, const int i) const {
+    assert(dir > 0 && dir < 4);
+    if (dir == X1DIR) return EdgeLength<X1DIR>(k, j, i);
+    else if (dir == X2DIR) return EdgeLength<X2DIR>(k, j, i);
+    return EdgeLength<X3DIR>(k, j, i);
+  }
+
+  //----------------------------------------
+  // FaceArea: Area of cell areas
+  //----------------------------------------
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const int k, const int j, const int i) const {
+    static_assert(dir > 0 && dir < 4);
+    if constexpr (dir == X1DIR) {
+      Real dOmega = Dx<X3DIR>() * (std::cos(Xf<X2DIR>(k, j, i)) - std::cos(Xf<X2DIR>(k, j+1, i)));
+      Real r = Xf<X1DIR>(k, j, i);
+      return r * r * dOmega;
+    } else if constexpr (dir == X2DIR) {
+      Real r0 = Xf<X1DIR>(k, j, i);
+      Real r1 = Xf<X1DIR>(k, j, i+1);
+      return (r1*r1*r1 - r0*r0*r0) * std::sin(Xf<X2DIR>(k, j, i)) * Dx<X3DIR>() / 3.0;
+    }
+    Real r0 = Xf<X1DIR>(k, j, i);
+    Real r1 = Xf<X1DIR>(k, j, i+1);
+    Real dcth = std::cos(Xf<X2DIR>(k, j, i)) - std::cos(Xf<X2DIR>(k, j+1, i));
+    return (r1*r1*r1 - r0*r0*r0) * dcth / 3.0;
+  }
+
+  //----------------------------------------
+  // CellVolume
+  //----------------------------------------
+  KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, const int i) const {
+    return FaceArea<X3DIR>(k, j, i) * Dx<X3DIR>();
+  }
+
+  KOKKOS_FORCEINLINE_FUNCTION
+  Real Volume(CellLevel cl, TopologicalElement el, const int k, const int j, const int i) {
+    using TE = TopologicalElement;
+    if (cl == CellLevel::same) {
+      if (el == TE::CC) return CellVolume(k, j, i);
+      else if (el == TE::F1) return FaceArea<X1DIR>(k, j, i);
+      else if (el == TE::F2) return FaceArea<X2DIR>(k, j, i);
+      else if (el == TE::F3) return FaceArea<X3DIR>(k, j, i);
+      else if (el == TE::E1) return EdgeLength<X1DIR>(k, j, i);
+      else if (el == TE::E2) return EdgeLength<X2DIR>(k, j, i);
+      else if (el == TE::E3) return EdgeLength<X3DIR>(k, j, i);
+      else if (el == TE::NN) return 1.0;
+    } else {
+      PARTHENON_FAIL("Have not yet implemented fine fields for UniformSpherical coordinates.");
+    }
+    PARTHENON_FAIL("If you reach this point, someone has added a new value to the the "
+                   "TopologicalElement enum.");
+    return 0.0;
+  }
+
+ private:
+};
+
+} // namespace parthenon
+
+#endif // COORDINATES_UNIFORM_SPHERICAL_HPP_
diff --git a/src/interface/swarm_device_context.hpp b/src/interface/swarm_device_context.hpp
index ed958126f36b..66a33aab6c0f 100644
--- a/src/interface/swarm_device_context.hpp
+++ b/src/interface/swarm_device_context.hpp
@@ -92,14 +92,14 @@ class SwarmDeviceContext {
   KOKKOS_INLINE_FUNCTION
   void Xtoijk(const Real &x, const Real &y, const Real &z, int &i, int &j, int &k) const {
     i = static_cast<int>(
-            std::floor((x - x_min_) / coords_.Dxc<CoordinateDirection::X1DIR>())) +
+            std::floor((x - x_min_) / coords_.Dx<CoordinateDirection::X1DIR>())) +
         ib_s_;
     j = (ndim_ > 1) ? static_cast<int>(std::floor(
-                          (y - y_min_) / coords_.Dxc<CoordinateDirection::X2DIR>())) +
+                          (y - y_min_) / coords_.Dx<CoordinateDirection::X2DIR>())) +
                           jb_s_
                     : jb_s_;
     k = (ndim_ > 2) ? static_cast<int>(std::floor(
-                          (z - z_min_) / coords_.Dxc<CoordinateDirection::X3DIR>())) +
+                          (z - z_min_) / coords_.Dx<CoordinateDirection::X3DIR>())) +
                           kb_s_
                     : kb_s_;
   }
diff --git a/src/prolong_restrict/pr_ops.hpp b/src/prolong_restrict/pr_ops.hpp
index ccf1e025c521..bc87f33f630a 100644
--- a/src/prolong_restrict/pr_ops.hpp
+++ b/src/prolong_restrict/pr_ops.hpp
@@ -443,9 +443,9 @@ struct ProlongateInternalTothAndRoe {
       const int dir1 = element_idx + 1;
       const int dir2 = (element_idx + 1) % 3 + 1;
       const int dir3 = (element_idx + 2) % 3 + 1;
-      const auto dx2 = std::pow(coarse_coords.DxcFA(dir1, k, j, i), 2);
-      const auto dy2 = std::pow(coarse_coords.DxcFA(dir2, k, j, i), 2);
-      const auto dz2 = std::pow(coarse_coords.DxcFA(dir3, k, j, i), 2);
+      const auto dx2 = std::pow(coarse_coords.Dxc(dir1, k, j, i), 2);
+      const auto dy2 = std::pow(coarse_coords.Dxc(dir2, k, j, i), 2);
+      const auto dz2 = std::pow(coarse_coords.Dxc(dir3, k, j, i), 2);
       Vxyz *= 0.125 * dz2 / (dx2 + dz2);
       Wxyz *= 0.125 * dy2 / (dx2 + dy2);
 
diff --git a/src/utils/index_split.hpp b/src/utils/index_split.hpp
index e9a5f4b441bd..4767bcff8aee 100644
--- a/src/utils/index_split.hpp
+++ b/src/utils/index_split.hpp
@@ -65,6 +65,16 @@ class IndexSplit {
   IndexRange GetInnerBounds(const IndexRange &jb, const IndexRange &ib) const {
     return {ib.s, (ibe_entire_ + 1) * (jb.e - jb.s + 1) - (ibe_entire_ - ib.e) - 1};
   }
+
+  KOKKOS_FORCEINLINE_FUNCTION
+  int get_i(const int idx) const {
+    return idx % (ibe_entire_ + 1);
+  }
+  KOKKOS_FORCEINLINE_FUNCTION
+  int get_deltaj(const int idx) const {
+    return idx / (ibe_entire_ + 1);
+  }
+
   KOKKOS_INLINE_FUNCTION
   bool is_i_ghost(const int idx) const {
     const int ni = ibe_entire_ + 1;
diff --git a/tst/unit/CMakeLists.txt b/tst/unit/CMakeLists.txt
index b44e6a7986da..595c62447a00 100644
--- a/tst/unit/CMakeLists.txt
+++ b/tst/unit/CMakeLists.txt
@@ -18,6 +18,7 @@
 
 list(APPEND unit_tests_SOURCES
     test_concepts_lite.cpp
+    test_coordinates.cpp
     test_data_collection.cpp
     test_taskid.cpp
     test_tasklist.cpp
diff --git a/tst/unit/test_coordinates.cpp b/tst/unit/test_coordinates.cpp
new file mode 100644
index 000000000000..c1e35fbb7c9f
--- /dev/null
+++ b/tst/unit/test_coordinates.cpp
@@ -0,0 +1,94 @@
+//========================================================================================
+// (C) (or copyright) 2024. 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.
+//========================================================================================
+
+#include <array>
+#include <iostream>
+
+#include "coordinates/uniform_cartesian.hpp"
+#include "coordinates/uniform_spherical.hpp"
+
+#include "basic_types.hpp"
+#include "defs.hpp"
+#include "globals.hpp"
+#include "parameter_input.hpp"
+using Real = double;
+using parthenon::ParameterInput;
+using parthenon::RegionSize;
+using parthenon::X1DIR;
+using parthenon::X2DIR;
+using parthenon::X3DIR;
+using parthenon::UniformCartesian;
+using parthenon::UniformSpherical;
+
+#include <catch2/catch.hpp>
+
+int nghost_save = 888;
+
+TEST_CASE("Checking UniformCartesian") {
+  std::array<Real, 3> xrat{1.0, 1.0, 1.0};
+  ParameterInput pin;
+  nghost_save = parthenon::Globals::nghost;
+  parthenon::Globals::nghost = 2;
+  GIVEN("A coordinate object") {
+    std::array<Real, 3> xmin{0.1, -0.2, 0.3};
+    std::array<Real, 3> xmax{0.3, 0.1, 0.7};
+    std::array<int, 3> nx{10, 12, 14};
+    RegionSize rs(xmin, xmax, xrat, nx);
+    UniformCartesian c(rs, &pin);
+    REQUIRE(c.Dx<X1DIR>() == (xmax[0] - xmin[0]) / nx[0]);
+    REQUIRE(c.Dx<X2DIR>() == (xmax[1] - xmin[1]) / nx[1]);
+    REQUIRE(c.Dx<X3DIR>() == (xmax[2] - xmin[2]) / nx[2]);
+  }
+}
+
+TEST_CASE("Checking UniformSpherical") {
+  std::array<Real, 3> xrat{1.0, 1.0, 1.0};
+  ParameterInput pin;
+  parthenon::Globals::nghost = 2;
+  GIVEN("A coordinate object") {
+    const Real rout  = 3.5;
+    const Real rin = 0.0;
+    std::array<Real, 3> xmin{rin, 0.0, 0.0};
+    std::array<Real, 3> xmax{rout, M_PI, 2 * M_PI};
+    std::array<int, 3> nx{3, 5, 8};
+    RegionSize rs(xmin, xmax, xrat, nx);
+    UniformSpherical c(rs, &pin);
+    const auto istart = c.GetStartIndex();
+    Real dr = (xmax[0] - xmin[0]) / nx[0];
+    const auto cxmin = c.GetXmin();
+    std::cout << "cxmin = " << cxmin[0] << "  " << "nghost = " << parthenon::Globals::nghost << std::endl;
+    REQUIRE(std::abs(cxmin[0] - (xmin[0] - parthenon::Globals::nghost*dr)) < 1.e-14);
+    int i0 = 6 + istart[0];
+    Real r0 = xmin[0] + (i0 - istart[0]) * dr;
+    REQUIRE(c.Xf<X1DIR>(i0) == r0);
+
+    Real area = 0.0;
+    for (int j = istart[1]; j < istart[1] + nx[1]; j++) {
+      for (int k = istart[2]; k < istart[2] + nx[2]; k++) {
+        area += c.FaceArea<X1DIR>(k, j, i0);
+      }
+    }
+    REQUIRE(std::abs(area - 4.0 * M_PI * r0 * r0) / area < 1.e-13);
+
+    Real volume = 0.0;
+    for (int k = istart[2]; k < istart[2] + nx[2]; k++) {
+      for (int j = istart[1]; j < istart[1] + nx[1]; j++) {
+        for (int i = istart[0]; i < istart[0] + nx[0]; i++) {
+          volume += c.CellVolume(k, j, i);
+        }
+      }
+    }
+    REQUIRE(std::abs(volume - (4.0 / 3.0 * M_PI * (std::pow(rout,3) - std::pow(rin,3)))) / volume < 1.e-13);
+  }
+  parthenon::Globals::nghost = nghost_save;
+}
\ No newline at end of file