diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp
index 6149e7f988..d9900c479d 100644
--- a/Source/driver/Castro.cpp
+++ b/Source/driver/Castro.cpp
@@ -3639,21 +3639,10 @@ Castro::apply_tagging_restrictions(TagBoxArray& tags, [[maybe_unused]] Real time
         {
             const Real* problo = geomdata.ProbLo();
             const Real* probhi = geomdata.ProbHi();
-            const Real* dx = geomdata.CellSize();
 
-            Real loc[3] = {0.0};
-
-            loc[0] = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0];
-#if AMREX_SPACEDIM >= 2
-            loc[1] = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1];
-#endif
-#if AMREX_SPACEDIM == 3
-            loc[2] = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2];
-#endif
-
-            Real r = std::sqrt((loc[0] - problem::center[0]) * (loc[0] - problem::center[0]) +
-                               (loc[1] - problem::center[1]) * (loc[1] - problem::center[1]) +
-                               (loc[2] - problem::center[2]) * (loc[2] - problem::center[2]));
+            GpuArray<Real, 3> loc;
+            position(i, j, k, geomdata, loc);
+            Real r = distance(geomdata, loc);
 
             Real max_dist_lo = 0.0;
             Real max_dist_hi = 0.0;
@@ -4357,9 +4346,12 @@ Castro::define_new_center(const MultiFab& S, Real time)
     // Now broadcast to everyone else.
     ParallelDescriptor::Bcast(&problem::center[0], AMREX_SPACEDIM, owner);
 
-    // Make sure if R-Z that center stays exactly on axis
+    // Make sure if R-Z and SPHERICAL that center stays exactly on axis
     if ( Geom().IsRZ() ) {
       problem::center[0] = 0;
+    } else if ( Geom().IsSPHERICAL() ) {
+      problem::center[0] = 0;
+      problem::center[1] = 0;
     }
 
 }
diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H
index 016215aa32..a9f3a87fae 100644
--- a/Source/driver/Castro_util.H
+++ b/Source/driver/Castro_util.H
@@ -94,7 +94,6 @@ bool mom_flux_has_p (const int mom_dir, const int flux_dir, const int coord)
 }
 
 
-
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void position(int i, int j, int k,
               GeometryData const& geomdata, GpuArray<Real, 3>& loc,
@@ -147,7 +146,7 @@ void position(int i, int j, int k,
     }
 
     for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
-        loc[dir] = offset[dir] + idx[dir] * dx[dir];
+        loc[dir] = offset[dir] + idx[dir] * dx[dir] - problem::center[dir];
     }
 
     for (int dir = AMREX_SPACEDIM; dir < 3; ++dir) {
@@ -156,6 +155,28 @@ void position(int i, int j, int k,
 
 }
 
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+Real distance(GeometryData const& geomdata, GpuArray<Real, 3>& loc)
+{
+    // Returns the distance from the center provided with loc, the position.
+    // loc is the form of [x,y,z] in Cartesian, [r, z, phi] in cylindrical
+    // and [r, theta, phi] in spherical
+
+    const auto coord = geomdata.Coord();
+
+    if (coord == CoordSys::cartesian) {
+        return std::sqrt(loc[0]*loc[0] + loc[1]*loc[1] + loc[2]*loc[2]);
+    }
+
+    if (coord == CoordSys::RZ) {
+        return std::sqrt(loc[0]*loc[0] + loc[1]*loc[1]);
+    }
+
+    return loc[0];
+}
+
+
 namespace geometry_util
 {
 
diff --git a/Source/driver/Derive.H b/Source/driver/Derive.H
index 7dd302da16..0533c2d824 100644
--- a/Source/driver/Derive.H
+++ b/Source/driver/Derive.H
@@ -12,125 +12,125 @@ extern "C"
 
   void ca_derpres
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_dereint1
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_dereint2
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_derlogden
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_deruplusc
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_deruminusc
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_dersoundspeed
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_dergamma1
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_dermachnumber
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_derentropy
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
 #ifdef DIFFUSION
   void ca_dercond
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_derdiffcoeff
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_derdiffterm
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 #endif
 
   void ca_derenuc
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_derenuctimescale
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_dervel
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_dermagvel
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_dermaggrav
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_derradialvel
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_dercircvel
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_dermagmom
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_derangmomx (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
-                      const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+                      const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
                       amrex::Real time, const int* bcrec, int level);
 
   void ca_derangmomy (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
-                      const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+                      const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
                       amrex::Real time, const int* bcrec, int level);
 
   void ca_derangmomz (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
-                      const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+                      const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
                       amrex::Real time, const int* bcrec, int level);
 
   void ca_derkineng (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
-                     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+                     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
                      amrex::Real time, const int* bcrec, int level);
 
   void ca_dernull
@@ -144,53 +144,53 @@ extern "C"
 
   void ca_derspec
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_derabar
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_derye
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_dermagvort
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_derdivu
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_derstate
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
 #ifdef MHD
   void ca_dermagcenx
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_dermagceny
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_dermagcenz
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
   void ca_derdivb
     (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
-     const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+     const amrex::FArrayBox& datfab, const amrex::Geometry& geom,
      amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
 
 #endif
diff --git a/Source/driver/Derive.cpp b/Source/driver/Derive.cpp
index 3064fc3dd2..f1f5f70f2c 100644
--- a/Source/driver/Derive.cpp
+++ b/Source/driver/Derive.cpp
@@ -21,7 +21,7 @@ extern "C"
     // need to explicitly synchronize after GPU kernels.
 
     void ca_derpres(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                    const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                    const FArrayBox& datfab, const Geometry& /*geom*/,
                     Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -54,7 +54,7 @@ extern "C"
     }
 
     void ca_dereint1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                     const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                     const FArrayBox& datfab, const Geometry& /*geom*/,
                      Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -76,7 +76,7 @@ extern "C"
     }
 
     void ca_dereint2(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                     const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                     const FArrayBox& datfab, const Geometry& /*geom*/,
                      Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -92,7 +92,7 @@ extern "C"
     }
 
     void ca_derlogden(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                      const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                      const FArrayBox& datfab, const Geometry& /*geom*/,
                       Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -107,7 +107,7 @@ extern "C"
     }
 
     void ca_deruplusc(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                      const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                      const FArrayBox& datfab, const Geometry& /*geom*/,
                       Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -142,7 +142,7 @@ extern "C"
     }
 
     void ca_deruminusc(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                       const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                       const FArrayBox& datfab, const Geometry& /*geom*/,
                        Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -177,7 +177,7 @@ extern "C"
     }
 
     void ca_dersoundspeed(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                          const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                          const FArrayBox& datfab, const Geometry& /*geom*/,
                           Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -213,7 +213,7 @@ extern "C"
 
 
     void ca_dergamma1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                      const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                      const FArrayBox& datfab, const Geometry& /*geom*/,
                       Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -248,7 +248,7 @@ extern "C"
     }
 
     void ca_dermachnumber(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                          const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                          const FArrayBox& datfab, const Geometry& /*geom*/,
                           Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -286,7 +286,7 @@ extern "C"
     }
 
     void ca_derentropy(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                       const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                       const FArrayBox& datfab, const Geometry& /*geom*/,
                        Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -321,7 +321,7 @@ extern "C"
 
 #ifdef DIFFUSION
     void ca_dercond(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                    const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                    const FArrayBox& datfab, const Geometry& /*geom*/,
                     Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -333,7 +333,7 @@ extern "C"
     }
 
     void ca_derdiffcoeff(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                    const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                    const FArrayBox& datfab, const Geometry& /*geom*/,
                     Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -345,7 +345,7 @@ extern "C"
     }
 
     void ca_derdiffterm(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                        const FArrayBox& datfab, const Geometry& geomdata,
+                        const FArrayBox& datfab, const Geometry& geom,
                         Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -360,9 +360,9 @@ extern "C"
 
       fill_temp_cond(obx, dat, coeff_arr);
 
-      auto dx = geomdata.CellSizeArray();
-      auto problo = geomdata.ProbLoArray();
-      const int coord_type = geomdata.Coord();
+      auto dx = geom.CellSizeArray();
+      auto problo = geom.ProbLoArray();
+      const int coord_type = geom.Coord();
 
       amrex::ParallelFor(bx,
       [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -432,14 +432,14 @@ extern "C"
 
 #ifdef REACTIONS
     void ca_derenuctimescale(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                             const FArrayBox& datfab, const Geometry& geomdata,
+                             const FArrayBox& datfab, const Geometry& geom,
                              Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
       auto const dat = datfab.array();
       auto const der = derfab.array();
 
-      auto dx = geomdata.CellSizeArray();
+      auto dx = geom.CellSizeArray();
 
       Real dd = 0.0_rt;
 #if AMREX_SPACEDIM == 1
@@ -494,7 +494,7 @@ extern "C"
     }
 
     void ca_derenuc(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                    const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                    const FArrayBox& datfab, const Geometry& /*geom*/,
                     Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
       auto const dat = datfab.array();
@@ -512,7 +512,7 @@ extern "C"
 #endif
 
     void ca_dervel(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                   const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                   const FArrayBox& datfab, const Geometry& /*geom*/,
                    Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -528,7 +528,7 @@ extern "C"
 
 
     void ca_dermagvel(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                      const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                      const FArrayBox& datfab, const Geometry& /*geom*/,
                       Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -549,7 +549,7 @@ extern "C"
 
 
     void ca_dermaggrav(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                       const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                       const FArrayBox& datfab, const Geometry& /*geom*/,
                        Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -568,7 +568,7 @@ extern "C"
     }
 
     void ca_derradialvel(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                         const FArrayBox& datfab, const Geometry& geomdata,
+                         const FArrayBox& datfab, const Geometry& geom,
                          Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -577,24 +577,13 @@ extern "C"
       auto const dat = datfab.array();
       auto const der = derfab.array();
 
-      auto dx = geomdata.CellSizeArray();
-      auto problo = geomdata.ProbLoArray();
+      auto geomdata = geom.data();
 
       amrex::ParallelFor(bx,
       [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
       {
-
-        Real x = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0] - problem::center[0];
-#if AMREX_SPACEDIM >= 2
-        Real y = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1] - problem::center[1];
-#else
-        Real y = 0.0_rt;
-#endif
-#if AMREX_SPACEDIM == 3
-        Real z = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
-#else
-        Real z = 0.0_rt;
-#endif
+        GpuArray<Real, 3> loc;
+        position(i, j, k, geomdata, loc);
 
         if (domain_is_plane_parallel) {
 #if AMREX_SPACEDIM == 2
@@ -607,15 +596,15 @@ extern "C"
           // where e_r and e_phi are the cylindrical unit vectors
 
           // we need the distance in the x-y plane from the origin
-          Real r = std::sqrt(x*x + y*y);
-          der(i,j,k,0) = (dat(i,j,k,1)*x + dat(i,j,k,2)*y) / (dat(i,j,k,0)*r);
+          Real r = std::sqrt(loc[0]*loc[0] + loc[1]*loc[1]);
+          der(i,j,k,0) = (dat(i,j,k,1)*loc[0] + dat(i,j,k,2)*loc[1]) / (dat(i,j,k,0)*r);
 #endif
         } else {
-          Real r = std::sqrt(x*x + y*y + z*z);
+          Real r = distance(geomdata, loc);
 
-          der(i,j,k,0) = (dat(i,j,k,1)*x +
-                          dat(i,j,k,2)*y +
-                          dat(i,j,k,3)*z) / ( dat(i,j,k,0)*r );
+          der(i,j,k,0) = (dat(i,j,k,1)*loc[0] +
+                          dat(i,j,k,2)*loc[1] +
+                          dat(i,j,k,3)*loc[2]) / ( dat(i,j,k,0)*r );
         }
 
       });
@@ -623,7 +612,7 @@ extern "C"
 
 
     void ca_dercircvel(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                       const FArrayBox& datfab, const Geometry& geomdata,
+                       const FArrayBox& datfab, const Geometry& geom,
                        Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -632,24 +621,13 @@ extern "C"
       auto const dat = datfab.array();
       auto const der = derfab.array();
 
-      auto dx = geomdata.CellSizeArray();
-      auto problo = geomdata.ProbLoArray();
+      auto geomdata = geom.data();
 
       amrex::ParallelFor(bx,
       [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
       {
-
-        Real x = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0] - problem::center[0];
-#if AMREX_SPACEDIM >= 2
-        Real y = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1] - problem::center[1];
-#else
-        Real y = 0.0_rt;
-#endif
-#if AMREX_SPACEDIM == 3
-        Real z = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
-#else
-        Real z = 0.0_rt;
-#endif
+        GpuArray<Real, 3> loc;
+        position(i, j, k, geomdata, loc);
 
         if (domain_is_plane_parallel) {
 #if AMREX_SPACEDIM == 2
@@ -662,11 +640,11 @@ extern "C"
           // where e_r and e_phi are the cylindrical unit vectors
 
           // we need the distance in the x-y plane from the origin
-          Real r = std::sqrt(x*x + y*y);
-          der(i,j,k,0) = (-dat(i,j,k,1)*y + dat(i,j,k,2)*x) / (dat(i,j,k,0)*r);
+          Real r = std::sqrt(loc[0]*loc[0] + loc[1]*loc[1]);
+          der(i,j,k,0) = (-dat(i,j,k,1)*loc[1] + dat(i,j,k,2)*loc[0]) / (dat(i,j,k,0)*r);
 #endif
         } else {
-          Real r = std::sqrt(x*x + y*y + z*z);
+          Real r = distance(geomdata, loc);
 
           // we really mean just the velocity component that is
           // perpendicular to radial, and in general 3-d (e.g. a
@@ -676,9 +654,9 @@ extern "C"
                         dat(i,j,k,2)*dat(i,j,k,2) +
                         dat(i,j,k,3)*dat(i,j,k,3))/(dat(i,j,k,0)*dat(i,j,k,0));
 
-          Real vr = (dat(i,j,k,1)*x +
-                     dat(i,j,k,2)*y +
-                     dat(i,j,k,3)*z) / ( dat(i,j,k,0)*r );
+          Real vr = (dat(i,j,k,1)*loc[0] +
+                     dat(i,j,k,2)*loc[1] +
+                     dat(i,j,k,3)*loc[2]) / ( dat(i,j,k,0)*r );
 
           der(i,j,k,0) = std::sqrt(amrex::max(vtot2 - vr*vr, 0.0_rt));
         }
@@ -688,7 +666,7 @@ extern "C"
 
 
   void ca_dermagmom(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                    const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                    const FArrayBox& datfab, const Geometry& /*geom*/,
                     Real /*time*/, const int* /*bcrec*/, int /*level*/)
   {
 
@@ -707,13 +685,12 @@ extern "C"
   }
 
   void ca_derangmomx (const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                      const FArrayBox& datfab, const Geometry& geomdata,
+                      const FArrayBox& datfab, const Geometry& geom,
                       Real /*time*/, const int* /*bcrec*/, int /*level*/)
   {
 
     int idir = 0;
-    auto dx     = geomdata.CellSizeArray();
-    auto problo = geomdata.ProbLoArray();
+    auto geomdata = geom.data();
 
     auto const dat = datfab.array();
     auto const L = derfab.array();
@@ -721,27 +698,8 @@ extern "C"
     amrex::ParallelFor(bx,
                        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                        {
-                         Real loc[3];
-
-                         //loc calculated like sum_utils.cpp
-                         //This might be equivalent and more modular: position(i, j, k, geomdata, loc);
-                         loc[0] = problo[0] + (0.5_rt + i) * dx[0];
-
-#if AMREX_SPACEDIM >= 2
-                         loc[1] = problo[1] + (0.5_rt + j) * dx[1];
-#else
-                         loc[1] = 0.0_rt;
-#endif
-
-#if AMREX_SPACEDIM == 3
-                         loc[2] = problo[2] + (0.5_rt + k) * dx[2];
-#else
-                         loc[2] = 0.0_rt;
-#endif
-
-                         for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
-                           loc[dir] -= problem::center[dir];
-                         }
+                         GpuArray<Real, 3> loc;
+                         position(i, j, k, geomdata, loc);
 
                          // Explicitly computing only the required cross-product as in inertial_to_rotational_velocity
                          if (idir == 0) { // cross_product(loc, mom): ang_mom(1)->x)
@@ -759,13 +717,12 @@ extern "C"
   }
 
   void ca_derangmomy (const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                      const FArrayBox& datfab, const Geometry& geomdata,
+                      const FArrayBox& datfab, const Geometry& geom,
                       Real /*time*/, const int* /*bcrec*/, int /*level*/)
   {
 
     int idir = 1;
-    auto dx     = geomdata.CellSizeArray();
-    auto problo = geomdata.ProbLoArray();
+    auto geomdata = geom.data();
 
     auto const dat = datfab.array();
     auto const L = derfab.array();
@@ -773,24 +730,8 @@ extern "C"
     amrex::ParallelFor(bx,
                        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                        {
-                         Real loc[3];
-
-                         loc[0] = problo[0] + (0.5_rt + i) * dx[0];
-
-#if AMREX_SPACEDIM >= 2
-                         loc[1] = problo[1] + (0.5_rt + j) * dx[1];
-#else
-                         loc[1] = 0.0_rt;
-#endif
-
-#if AMREX_SPACEDIM == 3
-                         loc[2] = problo[2] + (0.5_rt + k) * dx[2];
-#else
-                         loc[2] = 0.0_rt;
-#endif
-                         for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
-                           loc[dir] -= problem::center[dir];
-                         }
+                         GpuArray<Real, 3> loc;
+                         position(i, j, k, geomdata, loc);
 
                          if (idir == 0) { // cross_product(loc, mom): ang_mom(1)->x)
                            L(i,j,k,0) = loc[1] * dat(i,j,k,3) - loc[2] * dat(i,j,k,2);
@@ -807,13 +748,12 @@ extern "C"
   }
 
   void ca_derangmomz (const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                      const FArrayBox& datfab, const Geometry& geomdata,
+                      const FArrayBox& datfab, const Geometry& geom,
                       Real /*time*/, const int* /*bcrec*/, int /*level*/)
   {
 
     int idir = 2;
-    auto dx     = geomdata.CellSizeArray();
-    auto problo = geomdata.ProbLoArray();
+    auto geomdata = geom.data();
 
     auto const dat = datfab.array();
     auto const L = derfab.array();
@@ -821,25 +761,8 @@ extern "C"
     amrex::ParallelFor(bx,
                        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                        {
-                         Real loc[3];
-
-                         loc[0] = problo[0] + (0.5_rt + i) * dx[0];
-
-#if AMREX_SPACEDIM >= 2
-                         loc[1] = problo[1] + (0.5_rt + j) * dx[1];
-#else
-                         loc[1] = 0.0_rt;
-#endif
-
-#if AMREX_SPACEDIM == 3
-                         loc[2] = problo[2] + (0.5_rt + k) * dx[2];
-#else
-                         loc[2] = 0.0_rt;
-#endif
-
-                         for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
-                           loc[dir] -= problem::center[dir];
-                         }
+                         GpuArray<Real, 3> loc;
+                         position(i, j, k, geomdata, loc);
 
                          if (idir == 0) { // cross_product(loc, mom): ang_mom(1)->x)
                            L(i,j,k,0) = loc[1] * dat(i,j,k,3) - loc[2] * dat(i,j,k,2);
@@ -856,7 +779,7 @@ extern "C"
   }
 
   void ca_derkineng (const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                     const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                     const FArrayBox& datfab, const Geometry& /*geom*/,
                      Real /*time*/, const int* /*bcrec*/, int /*level*/)
   {
 
@@ -887,7 +810,7 @@ extern "C"
   }
 
   void ca_derspec(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                  const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                  const FArrayBox& datfab, const Geometry& /*geom*/,
                   Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -903,7 +826,7 @@ extern "C"
 
 
   void ca_derye(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                const FArrayBox& datfab, const Geometry& /*geom*/,
                 Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -925,7 +848,7 @@ extern "C"
     }
 
   void ca_derabar(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                  const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                  const FArrayBox& datfab, const Geometry& /*geom*/,
                   Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
@@ -947,19 +870,19 @@ extern "C"
     }
 
   void ca_dermagvort(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                     const FArrayBox& datfab, const Geometry& geomdata,
+                     const FArrayBox& datfab, const Geometry& geom,
                      Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
       auto const dat = datfab.array();
       auto const der = derfab.array();
 
-      auto dx = geomdata.CellSizeArray();
+      auto dx = geom.CellSizeArray();
 
-      const int coord_type = geomdata.Coord();
+      const int coord_type = geom.Coord();
 
 #if AMREX_SPACEDIM == 2
-      auto problo = geomdata.ProbLoArray();
+      auto problo = geom.ProbLoArray();
 #endif
 
       amrex::ParallelFor(bx,
@@ -1046,18 +969,18 @@ extern "C"
     }
 
   void ca_derdivu(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                         const FArrayBox& datfab, const Geometry& geomdata,
+                         const FArrayBox& datfab, const Geometry& geom,
                          Real /*time*/, const int* /*bcrec*/, int /*level*/)
     {
 
       auto const dat = datfab.array();
       auto const der = derfab.array();
 
-      auto dx = geomdata.CellSizeArray();
+      auto dx = geom.CellSizeArray();
 
-      auto problo = geomdata.ProbLoArray();
+      auto problo = geom.ProbLoArray();
 
-      const int coord_type = geomdata.Coord();
+      const int coord_type = geom.Coord();
 
       amrex::ParallelFor(bx,
       [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -1112,7 +1035,7 @@ extern "C"
   }
 
   void ca_derstate(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                   const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                   const FArrayBox& datfab, const Geometry& /*geom*/,
                    Real /*time*/, const int* /*bcrec*/, int /*level*/)
   {
 
@@ -1138,7 +1061,7 @@ extern "C"
 
 #ifdef MHD
   void ca_dermagcenx(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                     const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                     const FArrayBox& datfab, const Geometry& /*geom*/,
                      Real /*time*/, const int* /*bcrec*/, int /*level*/)
   {
 
@@ -1154,7 +1077,7 @@ extern "C"
   }
 
   void ca_dermagceny(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                     const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                     const FArrayBox& datfab, const Geometry& /*geom*/,
                      Real /*time*/, const int* /*bcrec*/, int /*level*/)
   {
 
@@ -1170,7 +1093,7 @@ extern "C"
   }
 
   void ca_dermagcenz(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                     const FArrayBox& datfab, const Geometry& /*geomdata*/,
+                     const FArrayBox& datfab, const Geometry& /*geom*/,
                      Real /*time*/, const int* /*bcrec*/, int /*level*/)
   {
 
@@ -1186,14 +1109,14 @@ extern "C"
   }
 
   void ca_derdivb(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
-                  const FArrayBox& datfab, const Geometry& geomdata,
+                  const FArrayBox& datfab, const Geometry& geom,
                   Real /*time*/, const int* /*bcrec*/, int /*level*/)
   {
 
     auto const dat = datfab.array();
     auto const der = derfab.array();
 
-    auto dx = geomdata.CellSizeArray();
+    auto dx = geom.CellSizeArray();
 
     amrex::ParallelFor(bx,
     [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
diff --git a/Source/driver/sum_utils.cpp b/Source/driver/sum_utils.cpp
index 6423b50b70..82d59bcf9f 100644
--- a/Source/driver/sum_utils.cpp
+++ b/Source/driver/sum_utils.cpp
@@ -384,10 +384,6 @@ Castro::gwstrain (Real time,
             GpuArray<Real, 3> r;
             position(i, j, k, geomdata, r);
 
-            for (int n = 0; n < 3; ++n) {
-                r[n] -= problem::center[n];
-            }
-
             Real rhoInv;
             if (rho(i,j,k) * maskFactor > 0.0_rt) {
                 rhoInv = 1.0_rt / rho(i,j,k);
diff --git a/Source/gravity/Castro_gravity.cpp b/Source/gravity/Castro_gravity.cpp
index 395e0e675e..cbada962bb 100644
--- a/Source/gravity/Castro_gravity.cpp
+++ b/Source/gravity/Castro_gravity.cpp
@@ -375,10 +375,7 @@ void Castro::construct_old_gravity_source(MultiFab& source, MultiFab& state_in,
 
 #ifdef HYBRID_MOMENTUM
             GpuArray<Real, 3> loc;
-            for (int n = 0; n < 3; ++n) {
-                position(i, j, k, geomdata, loc);
-                loc[n] -= problem::center[n];
-            }
+            position(i, j, k, geomdata, loc);
 
             GpuArray<Real, 3> hybrid_src;
 
@@ -574,9 +571,6 @@ void Castro::construct_new_gravity_source(MultiFab& source, MultiFab& state_old,
 #ifdef HYBRID_MOMENTUM
                 GpuArray<Real, 3> loc;
                 position(i, j, k, geomdata, loc);
-                for (int n = 0; n < 3; ++n) {
-                    loc[n] -= problem::center[n];
-                }
 
                 GpuArray<Real, 3> hybrid_src;
 
diff --git a/Source/gravity/Gravity.cpp b/Source/gravity/Gravity.cpp
index 82003be86f..a9cc9dbab1 100644
--- a/Source/gravity/Gravity.cpp
+++ b/Source/gravity/Gravity.cpp
@@ -1333,7 +1333,7 @@ Gravity::interpolate_monopole_grav(int level, RealVector& radial_grav, MultiFab&
     const auto dx = geom.CellSizeArray();
     const Real dr = dx[0] / static_cast<Real>(gravity::drdxfac);
 
-    const auto problo = geom.ProbLoArray();
+    const auto geomdata = geom.data();
 
 #ifdef _OPENMP
 #pragma omp parallel
@@ -1352,22 +1352,8 @@ Gravity::interpolate_monopole_grav(int level, RealVector& radial_grav, MultiFab&
         [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             GpuArray<Real, 3> loc;
-
-            loc[0] = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0] - problem::center[0];
-
-#if AMREX_SPACEDIM >= 2
-            loc[1] = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1] - problem::center[1];
-#else
-            loc[1] = 0.0_rt;
-#endif
-
-#if AMREX_SPACEDIM == 3
-            loc[2] = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
-#else
-            loc[2] = 0.0_rt;
-#endif
-
-            Real r = std::sqrt(loc[0] * loc[0] + loc[1] * loc[1] + loc[2] * loc[2]);
+            position(i, j, k, geomdata, loc);
+            Real r = distance(geomdata, loc);
 
             int index = static_cast<int>(r / dr);
 
@@ -1457,6 +1443,7 @@ Gravity::compute_radial_mass(const Box& bx,
     Real drinv = 1.0_rt / dr;
 
     const int coord_type = geom.Coord();
+    const auto geomdata = geom.data();
 
     AMREX_ALWAYS_ASSERT(coord_type >= 0 && coord_type <= 2);
 
@@ -1497,16 +1484,14 @@ Gravity::compute_radial_mass(const Box& bx,
     amrex::ParallelFor(bx,
     [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
-        Real xc = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0] - problem::center[0];
         Real lo_i = problo[0] + static_cast<Real>(i) * dx[0] - problem::center[0];
-
-        Real yc = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1] - problem::center[1];
         Real lo_j = problo[1] + static_cast<Real>(j) * dx[1] - problem::center[1];
-
-        Real zc = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
         Real lo_k = problo[2] + static_cast<Real>(k) * dx[2] - problem::center[2];
 
-        Real r = std::sqrt(xc * xc + yc * yc + zc * zc);
+        GpuArray<Real ,3> loc;
+        position(i, j, k, geomdata, loc);
+        Real r = distance(geomdata, loc);
+
         int index = static_cast<int>(r * drinv);
 
         // We may be coming in here with a masked out zone (in a zone on a coarse
diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp
index d443bbe060..70b2d74008 100644
--- a/Source/hydro/Castro_ctu_hydro.cpp
+++ b/Source/hydro/Castro_ctu_hydro.cpp
@@ -1212,9 +1212,6 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)  // NOLINT(readability-co
 
           position(i, j, k, geomdata, loc);
 
-          for (int dir = 0; dir < AMREX_SPACEDIM; ++dir)
-              loc[dir] -= problem::center[dir];
-
           Real R = amrex::max(std::sqrt(loc[0] * loc[0] + loc[1] * loc[1]),
                               std::numeric_limits<Real>::min());
           Real RInv = 1.0_rt / R;
diff --git a/Source/hydro/Castro_hybrid.cpp b/Source/hydro/Castro_hybrid.cpp
index 0191587b6f..8dfea244c9 100644
--- a/Source/hydro/Castro_hybrid.cpp
+++ b/Source/hydro/Castro_hybrid.cpp
@@ -107,9 +107,6 @@ Castro::fill_hybrid_hydro_source(MultiFab& sources, const MultiFab& state_in, Re
 
             position(i, j, k, geomdata, loc);
 
-            loc[0] -= problem::center[0];
-            loc[1] -= problem::center[1];
-
             Real R = amrex::max(std::sqrt(loc[0] * loc[0] + loc[1] * loc[1]),
                                 std::numeric_limits<Real>::min());
 
@@ -150,21 +147,19 @@ Castro::linear_to_hybrid_momentum(MultiFab& state_in, int ng)
 
             position(i, j, k, geomdata, loc);
 
-            for (int dir = 0; dir < AMREX_SPACEDIM; ++dir)
-                loc[dir] -= problem::center[dir];
-
             GpuArray<Real, 3> linear_mom;
 
-            for (int dir = 0; dir < 3; ++dir)
+            for (int dir = 0; dir < 3; ++dir) {
                 linear_mom[dir] = u(i,j,k,UMX+dir);
+            }
 
             GpuArray<Real, 3> hybrid_mom;
 
             linear_to_hybrid(loc, linear_mom, hybrid_mom);
 
-            for (int dir = 0; dir < 3; ++dir)
+            for (int dir = 0; dir < 3; ++dir) {
                 u(i,j,k,UMR+dir) = hybrid_mom[dir];
-
+            }
         });
     }
 }
@@ -196,21 +191,19 @@ Castro::hybrid_to_linear_momentum(MultiFab& state_in, int ng)
 
             position(i, j, k, geomdata, loc);
 
-            for (int dir = 0; dir < AMREX_SPACEDIM; ++dir)
-                loc[dir] -= problem::center[dir];
-
             GpuArray<Real, 3> hybrid_mom;
 
-            for (int dir = 0; dir < 3; ++dir)
+            for (int dir = 0; dir < 3; ++dir) {
                 hybrid_mom[dir] = u(i,j,k,UMR+dir);
+            }
 
             GpuArray<Real, 3> linear_mom;
 
             hybrid_to_linear(loc, hybrid_mom, linear_mom);
 
-            for (int dir = 0; dir < 3; ++dir)
+            for (int dir = 0; dir < 3; ++dir) {
                 u(i,j,k,UMX+dir) = linear_mom[dir];
-
+            }
         });
     }
 }
diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp
index 75899d67a1..d4f8747a07 100644
--- a/Source/hydro/Castro_mol_hydro.cpp
+++ b/Source/hydro/Castro_mol_hydro.cpp
@@ -636,9 +636,6 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update)
 
           position(i, j, k, geomdata, loc);
 
-          for (int dir = 0; dir < AMREX_SPACEDIM; ++dir)
-            loc[dir] -= problem::center[dir];
-
           Real R = amrex::max(std::sqrt(loc[0] * loc[0] + loc[1] * loc[1]),
                               std::numeric_limits<Real>::min());
           Real RInv = 1.0_rt / R;
diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp
index 515778c4d1..f87386d6dd 100644
--- a/Source/hydro/advection_util.cpp
+++ b/Source/hydro/advection_util.cpp
@@ -726,10 +726,6 @@ Castro::do_enforce_minimum_density(const Box& bx,
 
       position(i, j, k, geomdata, loc);
 
-      for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
-        loc[dir] -= problem::center[dir];
-      }
-
       GpuArray<Real, 3> linear_mom;
 
       for (int dir = 0; dir < 3; ++dir) {
diff --git a/Source/hydro/hybrid.H b/Source/hydro/hybrid.H
index 0188a1536e..632947f61f 100644
--- a/Source/hydro/hybrid.H
+++ b/Source/hydro/hybrid.H
@@ -123,13 +123,11 @@ void compute_hybrid_flux(const GpuArray<Real, NGDNV>& state, const GeometryData&
 
     position(i, j, k, geomdata, loc, ccx, ccy, ccz);
 
-    for (int dir = 0; dir < AMREX_SPACEDIM; ++dir)
-        loc[dir] -= problem::center[dir];
-
     GpuArray<Real, 3> linear_mom;
 
-    for (int dir = 0; dir < 3; ++dir)
+    for (int dir = 0; dir < 3; ++dir) {
         linear_mom[dir] = state[GDRHO] * state[GDU + dir];
+    }
 
     GpuArray<Real, 3> hybrid_mom;
 
diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp
index b7a422f1dd..69625998ec 100644
--- a/Source/reactions/Castro_react.cpp
+++ b/Source/reactions/Castro_react.cpp
@@ -208,7 +208,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra
 
         const auto dx = geom.CellSizeArray();
 #ifdef MODEL_PARSER
-        const auto problo = geom.ProbLoArray();
+        const auto geomdata = geom.data();
 #endif
 
 #if defined(AMREX_USE_GPU)
@@ -270,22 +270,15 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra
 
 #ifdef MODEL_PARSER
             if (drive_initial_convection) {
-                Real rr[3] = {0.0_rt};
-
-                rr[0] = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt) - problem::center[0];
-#if AMREX_SPACEDIM >= 2
-                rr[1] = problo[1] + dx[1] * (static_cast<Real>(j) + 0.5_rt) - problem::center[1];
-#endif
-#if AMREX_SPACEDIM == 3
-                rr[2] = problo[2] + dx[2] * (static_cast<Real>(k) + 0.5_rt) - problem::center[2];
-#endif
+                GpuArray<Real, 3> rr;
+                position(i, j, k, geomdata, rr);
 
                 Real dist;
 
                 if (domain_is_plane_parallel) {
                     dist = rr[AMREX_SPACEDIM-1];
                 } else {
-                    dist = std::sqrt(rr[0] * rr[0] + rr[1] * rr[1] + rr[2] * rr[2]);
+                    dist = distance(geomdata, rr);
                 }
 
                 burn_state.T_fixed = interpolate(dist, model::itemp);
@@ -564,7 +557,7 @@ Castro::react_state(Real time, Real dt)
         int lsdc_iteration = sdc_iteration;
 
         const auto dx = geom.CellSizeArray();
-        const auto problo = geom.ProbLoArray();
+        const auto geomdata = geom.data();
 
 #if defined(AMREX_USE_GPU)
         ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
@@ -618,22 +611,15 @@ Castro::react_state(Real time, Real dt)
 
 #ifdef MODEL_PARSER
             if (drive_initial_convection) {
-                Real rr[3] = {0.0_rt};
-
-                rr[0] = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt) - problem::center[0];
-#if AMREX_SPACEDIM >= 2
-                rr[1] = problo[1] + dx[1] * (static_cast<Real>(j) + 0.5_rt) - problem::center[1];
-#endif
-#if AMREX_SPACEDIM == 3
-                rr[2] = problo[2] + dx[2] * (static_cast<Real>(k) + 0.5_rt) - problem::center[2];
-#endif
+                GpuArray<Real, 3> rr;
+                position(i, j, k, geomdata, rr);
 
                 Real dist;
 
                 if (domain_is_plane_parallel) {
                     dist = rr[AMREX_SPACEDIM-1];
                 } else {
-                    dist = std::sqrt(rr[0] * rr[0] + rr[1] * rr[1] + rr[2] * rr[2]);
+                    dist = distance(geomdata, rr);
                 }
 
                 burn_state.T_fixed = interpolate(dist, model::itemp);
diff --git a/Source/rotation/Rotation.H b/Source/rotation/Rotation.H
index 878ae4a6b0..ef6c3e4129 100644
--- a/Source/rotation/Rotation.H
+++ b/Source/rotation/Rotation.H
@@ -129,10 +129,6 @@ inertial_to_rotational_velocity (const int i, const int j, const int k,
 
   position(i, j, k, geomdata, loc);
 
-  for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
-    loc[dir] -= problem::center[dir];
-  }
-
   auto omega = get_omega();
 
   // do the cross product Omega x loc
@@ -159,10 +155,6 @@ rotational_to_inertial_velocity (const int i, const int j, const int k,
 
   position(i, j, k, geomdata, loc);
 
-  for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
-    loc[dir] -= problem::center[dir];
-  }
-
   auto omega = get_omega();
 
   // do the cross product Omega x loc
diff --git a/Source/rotation/Rotation.cpp b/Source/rotation/Rotation.cpp
index 9a5f858fec..610b7a0e46 100644
--- a/Source/rotation/Rotation.cpp
+++ b/Source/rotation/Rotation.cpp
@@ -21,27 +21,14 @@ Castro::fill_rotational_psi(const Box& bx,
   auto omega = get_omega();
   Real denom = omega[0] * omega[0] + omega[1] * omega[1] + omega[2] * omega[2];
 
-  auto problo = geom.ProbLoArray();
-
-  auto dx = geom.CellSizeArray();
+  auto geomdata = geom.data();
 
   amrex::ParallelFor(bx,
   [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
   {
 
     GpuArray<Real, 3> r;
-
-    r[0] = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt) - problem::center[0];
-#if AMREX_SPACEDIM >= 2
-    r[1] = problo[1] + dx[1] * (static_cast<Real>(j) + 0.5_rt) - problem::center[1];
-#else
-    r[1] = 0.0_rt;
-#endif
-#if AMREX_SPACEDIM == 3
-    r[2] = problo[2] + dx[2] * (static_cast<Real>(k) + 0.5_rt) - problem::center[2];
-#else
-    r[2] = 0.0_rt;
-#endif
+    position(i, j, k, geomdata, r);
 
     if (denom != 0.0) {
         psi(i,j,k) = rotational_potential(r) / denom;
diff --git a/Source/rotation/rotation_sources.cpp b/Source/rotation/rotation_sources.cpp
index 4c034810d8..e09554d903 100644
--- a/Source/rotation/rotation_sources.cpp
+++ b/Source/rotation/rotation_sources.cpp
@@ -27,10 +27,6 @@ Castro::rsrc(const Box& bx,
     GpuArray<Real, 3> loc;
     position(i, j, k, geomdata, loc);
 
-    for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
-      loc[dir] -= problem::center[dir];
-    }
-
     Real rho = uold(i,j,k,URHO);
     Real rhoInv = 1.0_rt / rho;
 
@@ -241,10 +237,6 @@ Castro::corrrsrc(const Box& bx,
     GpuArray<Real, 3> loc;
     position(i, j, k, geomdata, loc);
 
-    for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
-      loc[dir] -= problem::center[dir];
-    }
-
     Real rhoo = uold(i,j,k,URHO);
     Real rhooinv = 1.0_rt / uold(i,j,k,URHO);
 
@@ -448,10 +440,6 @@ Castro::corrrsrc(const Box& bx,
 
               position(ie, je, ke, geomdata, loc, ccx, ccy, ccz);
 
-              for (int n = 0; n < AMREX_SPACEDIM; ++n) {
-                  loc[n] -= problem::center[n];
-              }
-
               GpuArray<Real, 3> temp_vel{};
               Real temp_Sr[3];
 
@@ -487,4 +475,3 @@ Castro::corrrsrc(const Box& bx,
   });
 
 }
-
diff --git a/Source/scf/scf_relax.cpp b/Source/scf/scf_relax.cpp
index 22935bb4e9..674a509f89 100644
--- a/Source/scf/scf_relax.cpp
+++ b/Source/scf/scf_relax.cpp
@@ -243,19 +243,11 @@ Castro::do_hscf_solve()
                 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
                 {
                     const auto *dx = geomdata.CellSize();
-                    const auto *problo = geomdata.ProbLo();
 
                     // The below assumes we are rotating on the z-axis.
 
-                    Real r[3] = {0.0};
-
-                    r[0] = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0] - problem::center[0];
-#if AMREX_SPACEDIM >= 2
-                    r[1] = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1] - problem::center[1];
-#endif
-#if AMREX_SPACEDIM == 3
-                    r[2] = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
-#endif
+                    GpuArray<Real, 3>r = {0.0};
+                    position(i, j, k, geomdata, r);
 
                     // Do a trilinear interpolation to find the contribution from
                     // this grid point. Limit so that only the nearest zone centers
@@ -376,19 +368,11 @@ Castro::do_hscf_solve()
                 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
                 {
                     const auto *dx = geomdata.CellSize();
-                    const auto *problo = geomdata.ProbLo();
 
                     // The below assumes we are rotating on the z-axis.
 
                     GpuArray<Real, 3> r = {0.0};
-
-                    r[0] = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0] - problem::center[0];
-#if AMREX_SPACEDIM >= 2
-                    r[1] = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1] - problem::center[1];
-#endif
-#if AMREX_SPACEDIM == 3
-                    r[2] = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
-#endif
+                    position(i, j, k, geomdata, r);
 
                     // Do a trilinear interpolation to find the contribution from
                     // this grid point. Limit so that only the nearest zone centers
@@ -455,18 +439,8 @@ Castro::do_hscf_solve()
                     // enthalpy + gravitational potential + rotational potential = const
                     // We already have the constant, so our goal is to construct the enthalpy field.
 
-                    const auto *dx = geomdata.CellSize();
-                    const auto *problo = geomdata.ProbLo();
-
                     GpuArray<Real, 3> r = {0.0};
-
-                    r[0] = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0] - problem::center[0];
-#if AMREX_SPACEDIM >= 2
-                    r[1] = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1] - problem::center[1];
-#endif
-#if AMREX_SPACEDIM == 3
-                    r[2] = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
-#endif
+                    position(i, j, k, geomdata, r);
 
                     enthalpy_arr(i,j,k) = bernoulli - phi_arr(i,j,k) - rotational_potential(r);
                 });
@@ -639,17 +613,8 @@ Castro::do_hscf_solve()
                 {
                     Real dM = 0.0, dK = 0.0, dU = 0.0, dE = 0.0;
 
-                    const auto* problo = geomdata.ProbLo();
-
                     GpuArray<Real, 3> r = {0.0};
-
-                    r[0] = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0] - problem::center[0];
-#if AMREX_SPACEDIM >= 2
-                    r[1] = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1] - problem::center[1];
-#endif
-#if AMREX_SPACEDIM == 3
-                    r[2] = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
-#endif
+                    position(i, j, k, geomdata, r);
 
                     if (state_arr(i,j,k,URHO) > 0.0)
                     {
diff --git a/Source/sources/Castro_sponge.cpp b/Source/sources/Castro_sponge.cpp
index 2109d63374..67be2b39a5 100644
--- a/Source/sources/Castro_sponge.cpp
+++ b/Source/sources/Castro_sponge.cpp
@@ -79,8 +79,7 @@ Castro::apply_sponge(const Box& bx,
     alpha = 0.0_rt;
   }
 
-  auto dx = geom.CellSizeArray();
-  auto problo = geom.ProbLoArray();
+  auto geomdata = geom.data();
 
   amrex::ParallelFor(bx,
   [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -89,20 +88,7 @@ Castro::apply_sponge(const Box& bx,
     Real src[NSRC] = {0.0};
 
     GpuArray<Real, 3> r;
-
-    r[0] = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0] - problem::center[0];
-
-#if AMREX_SPACEDIM >= 2
-    r[1] = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1] - problem::center[1];
-#else
-    r[1] = 0.0_rt;
-#endif
-
-#if AMREX_SPACEDIM == 3
-    r[2] = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
-#else
-    r[2] = 0.0_rt;
-#endif
+    position(i, j, k, geomdata, r);
 
     Real rho = state_in(i,j,k,URHO);
     Real rhoInv = 1.0_rt / rho;
@@ -124,7 +110,7 @@ Castro::apply_sponge(const Box& bx,
     Real sponge_factor = 0.0_rt;
 
     if (sponge_lower_radius >= 0.0_rt && sponge_upper_radius > sponge_lower_radius) {
-      Real rad = std::sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
+      Real rad = distance(geomdata, r);
 
       if (rad < sponge_lower_radius) {
         sponge_factor = sponge_lower_factor;