diff --git a/Exec/hydro_tests/rotating_torus/problem_initialize_state_data.H b/Exec/hydro_tests/rotating_torus/problem_initialize_state_data.H
index 30b70388b5..53a4630134 100644
--- a/Exec/hydro_tests/rotating_torus/problem_initialize_state_data.H
+++ b/Exec/hydro_tests/rotating_torus/problem_initialize_state_data.H
@@ -17,7 +17,7 @@ void problem_initialize_state_data (int i, int j, int k,
     const Real* probhi = geomdata.ProbHi();
 
 #ifdef ROTATION
-    auto omega = get_omega_vec(j);
+    auto omega = get_omega_vec(geomdata, j);
 #else
     // Provide a dummy value so that we can compile without rotation.
     Real omega[3] = {0.0_rt, 0.0_rt, 2.0_rt * M_PI};
diff --git a/Exec/science/wdmerger/Prob.cpp b/Exec/science/wdmerger/Prob.cpp
index 92303f67a5..78c7104463 100644
--- a/Exec/science/wdmerger/Prob.cpp
+++ b/Exec/science/wdmerger/Prob.cpp
@@ -798,6 +798,7 @@ Castro::update_relaxation(Real time, Real dt) {
     get_lagrange_points(mass_P, mass_S, com_P, com_S, L1, L2, L3);
 
     const auto dx = geom.CellSizeArray();
+    const auto coord = geom.Coord();
     GeometryData geomdata = geom.data();
 
     MultiFab& phi_new = get_new_data(PhiGrav_Type);
@@ -833,8 +834,8 @@ Castro::update_relaxation(Real time, Real dt) {
 
                 // Compute the effective potential.
 
-                auto omega = get_omega_vec(j);
-                Real phiEff = phi(i,j,k) + rotational_potential(r, omega);
+                auto omega = get_omega_vec(geomdata, j);
+                Real phiEff = phi(i,j,k) + rotational_potential(r, omega, coord);
 
                 for (int n = 0; n < 3; ++n) {
                     r[n] -= L1[n];
@@ -902,8 +903,8 @@ Castro::update_relaxation(Real time, Real dt) {
                 GpuArray<Real, 3> r;
                 position(i, j, k, geomdata, r);
 
-                auto omega = get_omega_vec(j);
-                Real phiEff = phi(i,j,k) + rotational_potential(r, omega);
+                auto omega = get_omega_vec(geomdata, j);
+                Real phiEff = phi(i,j,k) + rotational_potential(r, omega, coord);
 
                 Real done = 0.0_rt;
 
diff --git a/Exec/science/wdmerger/Problem_Derive.cpp b/Exec/science/wdmerger/Problem_Derive.cpp
index ef8e60b7c5..776b7078c5 100644
--- a/Exec/science/wdmerger/Problem_Derive.cpp
+++ b/Exec/science/wdmerger/Problem_Derive.cpp
@@ -419,6 +419,8 @@ void ca_derphieff(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/
 
     const auto dx = geom.CellSizeArray();
     const auto problo = geom.ProbLoArray();
+    const auto coord = geom.Coord();
+    const auto geomdata = geom.data();
 
     amrex::ParallelFor(bx,
     [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
@@ -438,9 +440,9 @@ void ca_derphieff(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/
 #else
         loc[2] = 0.0_rt;
 #endif
-        auto omega = get_omega_vec(j);
+        auto omega = get_omega_vec(geomdata, j);
 
-        der(i,j,k,0) = dat(i,j,k,0) + rotational_potential(loc, omega);
+        der(i,j,k,0) = dat(i,j,k,0) + rotational_potential(loc, omega, coord);
     });
 }
 
@@ -459,6 +461,8 @@ void ca_derphieffpm_p(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nco
 
     const auto dx = geom.CellSizeArray();
     const auto problo = geom.ProbLoArray();
+    const auto coord = geom.Coord();
+    const auto geomdata = geom.data();
 
     amrex::ParallelFor(bx,
     [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
@@ -493,9 +497,9 @@ void ca_derphieffpm_p(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nco
             loc[iloc] -= problem::center[iloc];
         }
 
-        auto omega = get_omega_vec(j);
+        auto omega = get_omega_vec(geomdata, j);
 
-        der(i,j,k,0) = -C::Gconst * problem::mass_P / r + rotational_potential(loc, omega);
+        der(i,j,k,0) = -C::Gconst * problem::mass_P / r + rotational_potential(loc, omega, coord);
     });
 }
 
@@ -511,6 +515,8 @@ void ca_derphieffpm_s(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nco
 
     const auto dx = geom.CellSizeArray();
     const auto problo = geom.ProbLoArray();
+    const auto coord = geom.Coord();
+    const auto geomdata = geom.data();
 
     amrex::ParallelFor(bx,
     [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
@@ -545,9 +551,9 @@ void ca_derphieffpm_s(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nco
             loc[iloc] -= problem::center[iloc];
         }
 
-        auto omega = get_omega_vec(j);
+        auto omega = get_omega_vec(geomdata, j);
 
-        der(i,j,k,0) = -C::Gconst * problem::mass_S / r + rotational_potential(loc, omega);
+        der(i,j,k,0) = -C::Gconst * problem::mass_S / r + rotational_potential(loc, omega, coord);
     });
 }
 
@@ -577,6 +583,8 @@ void ca_derrhophiRot(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncom
 
     const auto dx = geom.CellSizeArray();
     const auto problo = geom.ProbLoArray();
+    const auto coord = geom.Coord();
+    const auto geomdata = geom.data();
 
     amrex::ParallelFor(bx,
     [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
@@ -597,9 +605,9 @@ void ca_derrhophiRot(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncom
         loc[2] = 0.0_rt;
 #endif
 
-        auto omega = get_omega_vec(j);
+        auto omega = get_omega_vec(geomdata, j);
 
-        der(i,j,k,0) = dat(i,j,k,0) * rotational_potential(loc, omega);
+        der(i,j,k,0) = dat(i,j,k,0) * rotational_potential(loc, omega, coord);
     });
 }
 
diff --git a/Exec/science/wdmerger/problem_initialize_state_data.H b/Exec/science/wdmerger/problem_initialize_state_data.H
index 839f4abc25..f9f52414f8 100644
--- a/Exec/science/wdmerger/problem_initialize_state_data.H
+++ b/Exec/science/wdmerger/problem_initialize_state_data.H
@@ -16,7 +16,7 @@ void problem_initialize_state_data (int i, int j, int k,
     // or secondary (in which case interpolate from the respective model)
     // or if we are in an ambient zone.
 
-    auto omega = get_omega_vec(j);
+    auto omega = get_omega_vec(geomdata, j);
 
     GpuArray<Real, 3> loc;
     position(i, j, k, geomdata, loc);
diff --git a/Exec/science/wdmerger/wdmerger_util.H b/Exec/science/wdmerger/wdmerger_util.H
index d19a4c04cf..9a8c4fe327 100644
--- a/Exec/science/wdmerger/wdmerger_util.H
+++ b/Exec/science/wdmerger/wdmerger_util.H
@@ -89,9 +89,9 @@ Real stellar_mask (int i, int j, int k,
         loc[iloc] -= problem::center[iloc];
     }
 
-    auto omega = get_omega_vec(j);
+    auto omega = get_omega_vec(geomdata, j);
 
-    Real phi_rot = rotational_potential(loc, omega);
+    Real phi_rot = rotational_potential(loc, omega, geomdata.Coord());
 
     Real phi_p = -C::Gconst * problem::mass_P / r_P + phi_rot;
     Real phi_s = -C::Gconst * problem::mass_S / r_S + phi_rot;
diff --git a/Source/driver/sum_utils.cpp b/Source/driver/sum_utils.cpp
index bbdec07227..f19a02f2d6 100644
--- a/Source/driver/sum_utils.cpp
+++ b/Source/driver/sum_utils.cpp
@@ -400,7 +400,7 @@ Castro::gwstrain (Real time,
 
             GpuArray<Real, 3> pos{r};
 #ifdef ROTATION
-            auto omega = get_omega_vec(j);
+            auto omega = get_omega_vec(geomdata, j);
             pos = inertial_rotation(r, omega, time);
 #endif
 
@@ -432,7 +432,6 @@ Castro::gwstrain (Real time,
 
             GpuArray<Real, 3> inertial_g{g};
 #ifdef ROTATION
-            auto omega = get_omega_vec(j);
             inertial_g = inertial_rotation(g, omega, time);
 #endif
 
diff --git a/Source/rotation/Rotation.H b/Source/rotation/Rotation.H
index 514978f084..06f2a04490 100644
--- a/Source/rotation/Rotation.H
+++ b/Source/rotation/Rotation.H
@@ -29,20 +29,21 @@ Real get_omega()
 /// Return the omega vector corresponding to the current rotational period
 ///
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-GpuArray<Real, 3> get_omega_vec(const int j)
+GpuArray<Real, 3> get_omega_vec(const GeometryData& geomdata,
+                                const int j)
 {
     GpuArray<Real, 3> omega_vec = {0.0_rt};
 
-    const auto coord = geom.Coord();
-    const auto problo = geom.ProbLoArray();
-    const auto dx = geom.CellSizeArray();
+    const auto coord = geomdata.Coord();
+    const auto problo = geomdata.ProbLo(1);
+    const auto dtheta = geomdata.CellSize(1);
 
     auto omega = get_omega();
 
     // If coord == 2, i.e. Spherical2D, then it is assumed that rotational axis is in z-axis
     // then we convert back to r and theta direction.
     if (coord == 2) {
-        Real theta = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1];
+        Real theta = problo + (static_cast<Real>(j) + 0.5_rt) * dtheta;
         omega_vec[0] = std::cos(theta) * omega;
         omega_vec[1] = -std::sin(theta) * omega;
     } else {
@@ -58,13 +59,15 @@ GpuArray<Real, 3> get_omega_vec(const int j)
 ///
 /// @param r         distance from origin of rotation vector
 /// @param v         velocity
+/// @param omega     angular velocity vector array
+/// @param coord     coordinate system
 /// @param coriolis  do we include the Coriolis force
 /// @param Sr        rotational acceleration
 ///
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void
 rotational_acceleration(const GpuArray<Real, 3>& r, const GpuArray<Real, 3>& v,
-                        const GpuArray<Real, 3>& omega,
+                        const GpuArray<Real, 3>& omega, const int coord,
                         const bool coriolis, Real* Sr) {
 
   // Given a position and velocity, calculate
@@ -90,7 +93,6 @@ rotational_acceleration(const GpuArray<Real, 3>& r, const GpuArray<Real, 3>& v,
 
   bool c2 = (castro::rotation_include_coriolis == 1 && coriolis) ? true : false;
 
-  const auto coord = geom.Coord();
   auto r_vec = r;
 
   GpuArray<Real, 3> omega_cross_v;
@@ -125,14 +127,13 @@ rotational_acceleration(const GpuArray<Real, 3>& r, const GpuArray<Real, 3>& v,
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 Real
-rotational_potential(const GpuArray<Real, 3>& r, const GpuArray<Real, 3>& omega) {
+rotational_potential(const GpuArray<Real, 3>& r, const GpuArray<Real, 3>& omega,
+                     const int coord) {
 
   // Construct rotational potential, phi_R = -1/2 | omega x r |**2
   //
 
   Real phi = 0.0_rt;
-
-  const auto coord = geom.Coord();
   auto r_vec = r;
 
   if (coord == 2) {
@@ -176,7 +177,12 @@ inertial_to_rotational_velocity (const int i, const int j, const int k,
     loc[dir] -= problem::center[dir];
   }
 
-  auto omega = get_omega_vec(j);
+  if (geomdata.Coord() == 2) {
+      loc[1] = 0.0_rt;
+      loc[2] = 0.0_rt;
+  }
+
+  auto omega = get_omega_vec(geomdata, j);
 
   // do the cross product Omega x loc
   v[0] += -(omega[1]*loc[2] - omega[2]*loc[1]);
@@ -206,7 +212,12 @@ rotational_to_inertial_velocity (const int i, const int j, const int k,
     loc[dir] -= problem::center[dir];
   }
 
-  auto omega = get_omega_vec(j);
+  if (geomdata.Coord() == 2) {
+      loc[1] = 0.0_rt;
+      loc[2] = 0.0_rt;
+  }
+
+  auto omega = get_omega_vec(geomdata, j);
 
   // do the cross product Omega x loc
   v[0] += (omega[1]*loc[2] - omega[2]*loc[1]);
diff --git a/Source/rotation/Rotation.cpp b/Source/rotation/Rotation.cpp
index 3a270b3a3e..a495f49d1e 100644
--- a/Source/rotation/Rotation.cpp
+++ b/Source/rotation/Rotation.cpp
@@ -18,12 +18,13 @@ Castro::fill_rotational_psi(const Box& bx,
   // routine uniquely determines the rotation law. For the other
   // rotation laws, we would simply divide by v_0^2 or j_0^2 instead.
 
-  auto omega = get_omega();
+  const auto omega = get_omega();
   Real denom = omega * omega;
 
-  auto problo = geom.ProbLoArray();
-
-  auto dx = geom.CellSizeArray();
+  const auto problo = geom.ProbLoArray();
+  const auto dx = geom.CellSizeArray();
+  const auto coord = geom.Coord();
+  const auto geomdata = geom.data();
 
   amrex::ParallelFor(bx,
   [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -44,8 +45,8 @@ Castro::fill_rotational_psi(const Box& bx,
 #endif
 
     if (denom != 0.0) {
-        auto omega = get_omega_vec(j);
-        psi(i,j,k) = rotational_potential(r, omega) / denom;
+        auto omega = get_omega_vec(geomdata, j);
+        psi(i,j,k) = rotational_potential(r, omega, coord) / denom;
     }
     else {
         psi(i,j,k) = 0.0;
diff --git a/Source/rotation/rotation_sources.cpp b/Source/rotation/rotation_sources.cpp
index 4a6060d8ba..d116da5321 100644
--- a/Source/rotation/rotation_sources.cpp
+++ b/Source/rotation/rotation_sources.cpp
@@ -12,6 +12,7 @@ Castro::rsrc(const Box& bx,
              const Real dt) {
 
   GeometryData geomdata = geom.data();
+  const auto coord = geomdata.Coord();
 
   amrex::ParallelFor(bx,
   [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -31,7 +32,7 @@ Castro::rsrc(const Box& bx,
       loc[dir] -= problem::center[dir];
     }
 
-    auto omega = get_omega_vec(j);
+    auto omega = get_omega_vec(geomdata, j);
 
     Real rho = uold(i,j,k,URHO);
     Real rhoInv = 1.0_rt / rho;
@@ -49,7 +50,7 @@ Castro::rsrc(const Box& bx,
     v[2] = uold(i,j,k,UMZ) * rhoInv;
 
     bool coriolis = true;
-    rotational_acceleration(loc, v, omega, coriolis, Sr);
+    rotational_acceleration(loc, v, omega, coord, coriolis, Sr);
 
     for (auto& e : Sr) {
         e *= rho;
@@ -169,6 +170,7 @@ Castro::corrrsrc(const Box& bx,
   // is the new time at time-level n+1.
 
   GeometryData geomdata = geom.data();
+  const auto coord = geomdata.Coord();
 
   Real hdtInv = 0.5_rt / dt;
 
@@ -200,7 +202,7 @@ Castro::corrrsrc(const Box& bx,
       loc[dir] -= problem::center[dir];
     }
 
-    auto omega = get_omega_vec(j);
+    auto omega = get_omega_vec(geomdata, j);
 
     Real rhoo = uold(i,j,k,URHO);
     Real rhooinv = 1.0_rt / uold(i,j,k,URHO);
@@ -224,7 +226,7 @@ Castro::corrrsrc(const Box& bx,
     vold[2] = uold(i,j,k,UMZ) * rhooinv;
 
     bool coriolis = true;
-    rotational_acceleration(loc, vold, omega, coriolis, Sr_old);
+    rotational_acceleration(loc, vold, omega, coord, coriolis, Sr_old);
 
     for (auto& e : Sr_old) {
         e *= rhoo;
@@ -241,7 +243,7 @@ Castro::corrrsrc(const Box& bx,
     vnew[1] = unew(i,j,k,UMY) * rhoninv;
     vnew[2] = unew(i,j,k,UMZ) * rhoninv;
 
-    rotational_acceleration(loc, vnew, omega, coriolis, Sr_new);
+    rotational_acceleration(loc, vnew, omega, coord, coriolis, Sr_new);
 
     for (auto& e : Sr_new) {
         e *= rhon;
@@ -266,7 +268,7 @@ Castro::corrrsrc(const Box& bx,
 
       Real acc[3];
       coriolis = false;
-      rotational_acceleration(loc, vnew, omega, coriolis, acc);
+      rotational_acceleration(loc, vnew, omega, coord, coriolis, acc);
 
       Real new_mom_tmp[3];
       for (int n = 0; n < 3; n++) {
@@ -358,7 +360,7 @@ Castro::corrrsrc(const Box& bx,
 
       Real acc[3];
       coriolis = true;
-      rotational_acceleration(loc, vnew, omega, coriolis, acc);
+      rotational_acceleration(loc, vnew, omega, coord, coriolis, acc);
 
       Sr_new[0] = rhon * acc[0];
       Sr_new[1] = rhon * acc[1];
@@ -416,7 +418,7 @@ Castro::corrrsrc(const Box& bx,
               Real temp_Sr[3];
 
               coriolis = false;
-              rotational_acceleration(loc, temp_vel, omega, coriolis, temp_Sr);
+              rotational_acceleration(loc, temp_vel, omega, coord, coriolis, temp_Sr);
 
               edge_Sr[dir][edge] = temp_Sr[dir];
 
diff --git a/Source/scf/scf_relax.cpp b/Source/scf/scf_relax.cpp
index 0eecfe47f9..edd19b3c9a 100644
--- a/Source/scf/scf_relax.cpp
+++ b/Source/scf/scf_relax.cpp
@@ -244,6 +244,7 @@ Castro::do_hscf_solve()
                 {
                     const auto *dx = geomdata.CellSize();
                     const auto *problo = geomdata.ProbLo();
+                    const auto coord = geomdata.Coord();
 
                     // The below assumes we are rotating on the z-axis.
 
@@ -413,8 +414,8 @@ Castro::do_hscf_solve()
                         scale = (1.0_rt - rr[0]) * (1.0_rt - rr[1]) * (1.0_rt - rr[2]);
                     }
 
-                    auto omega = get_omega_vec(j);
-                    Real bernoulli_zone = scale * (phi_arr(i,j,k) + rotational_potential(r, omega));
+                    auto omega = get_omega_vec(geomdata, j);
+                    Real bernoulli_zone = scale * (phi_arr(i,j,k) + rotational_potential(r, omega, coord));
 
                     return {bernoulli_zone};
                 });
@@ -458,6 +459,7 @@ Castro::do_hscf_solve()
 
                     const auto *dx = geomdata.CellSize();
                     const auto *problo = geomdata.ProbLo();
+                    const auto coord = geomdata.Coord();
 
                     GpuArray<Real, 3> r = {0.0};
 
@@ -469,8 +471,8 @@ Castro::do_hscf_solve()
                     r[2] = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
 #endif
 
-                    auto omega = get_omega_vec(j);
-                    enthalpy_arr(i,j,k) = bernoulli - phi_arr(i,j,k) - rotational_potential(r, omega);
+                    auto omega = get_omega_vec(geomdata, j);
+                    enthalpy_arr(i,j,k) = bernoulli - phi_arr(i,j,k) - rotational_potential(r, omega, coord);
                 });
 
             }
@@ -642,6 +644,7 @@ Castro::do_hscf_solve()
                     Real dM = 0.0, dK = 0.0, dU = 0.0, dE = 0.0;
 
                     const auto* problo = geomdata.ProbLo();
+                    const auto coord = geomdata.Coord();
 
                     GpuArray<Real, 3> r = {0.0};
 
@@ -657,8 +660,8 @@ Castro::do_hscf_solve()
                     {
                         dM = state_arr(i,j,k,URHO) * dV;
 
-                        auto omega = get_omega_vec(j);
-                        dK = rotational_potential(r, omega) * dM;
+                        auto omega = get_omega_vec(geomdata, j);
+                        dK = rotational_potential(r, omega, coord) * dM;
 
                         dU = 0.5_rt * phi_arr(i,j,k) * dM;