diff --git a/include/DirectSolver/DirectSolverGive/applySymmetryShift.inl b/include/DirectSolver/DirectSolverGive/applySymmetryShift.inl
index 959d1c27..56f48a9e 100644
--- a/include/DirectSolver/DirectSolverGive/applySymmetryShift.inl
+++ b/include/DirectSolver/DirectSolverGive/applySymmetryShift.inl
@@ -24,6 +24,7 @@ void DirectSolverGive<LevelCacheType>::applySymmetryShiftInnerBoundary(Vector<do
 
             for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) {
                 const double theta = grid.theta(i_theta);
+
                 /* -------------------------- */
                 /* Node on the inner boundary */
                 /* -------------------------- */
@@ -37,12 +38,15 @@ void DirectSolverGive<LevelCacheType>::applySymmetryShiftInnerBoundary(Vector<do
                 k1 = grid.angularSpacing(i_theta - 1);
                 k2 = grid.angularSpacing(i_theta);
 
+                const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
+                const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
+
                 coeff2 = 0.5 * (k1 + k2) / h2;
 
                 /* Fill x(i+1,j) */
                 x(grid.index(i_r + 1, i_theta)) -= -coeff2 * arr * x(grid.index(i_r, i_theta)) /* Left */
-                                                   + 0.25 * art * x(grid.index(i_r, i_theta + 1)) /* Top Left */
-                                                   - 0.25 * art * x(grid.index(i_r, i_theta - 1)); /* Bottom Left */
+                                                   + 0.25 * art * x(grid.index(i_r, i_theta_P1)) /* Top Left */
+                                                   - 0.25 * art * x(grid.index(i_r, i_theta_M1)); /* Bottom Left */
 
                 /* --------------------------- */
                 /* Node next to inner boundary */
@@ -62,9 +66,9 @@ void DirectSolverGive<LevelCacheType>::applySymmetryShiftInnerBoundary(Vector<do
                 /* Fill x(i,j) */
                 x(grid.index(i_r, i_theta)) -= -coeff1 * arr * x(grid.index(i_r - 1, i_theta)); /* Left */
                 /* Fill x(i,j-1) */
-                x(grid.index(i_r, i_theta - 1)) -= +0.25 * art * x(grid.index(i_r - 1, i_theta)); /* Top Left */
+                x(grid.index(i_r, i_theta_M1)) -= +0.25 * art * x(grid.index(i_r - 1, i_theta)); /* Top Left */
                 /* Fill x(i,j+1) */
-                x(grid.index(i_r, i_theta + 1)) -= -0.25 * art * x(grid.index(i_r - 1, i_theta)); /* Bottom Left */
+                x(grid.index(i_r, i_theta_P1)) -= -0.25 * art * x(grid.index(i_r - 1, i_theta)); /* Bottom Left */
             }
         });
 }
@@ -100,14 +104,17 @@ void DirectSolverGive<LevelCacheType>::applySymmetryShiftOuterBoundary(Vector<do
                 k1 = grid.angularSpacing(i_theta - 1);
                 k2 = grid.angularSpacing(i_theta);
 
+                const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
+                const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
+
                 coeff2 = 0.5 * (k1 + k2) / h2;
 
                 /* Fill result(i,j) */
                 x(grid.index(i_r, i_theta)) -= -coeff2 * arr * x(grid.index(i_r + 1, i_theta)); /* Right */
                 /* Fill result(i,j-1) */
-                x(grid.index(i_r, i_theta - 1)) -= -0.25 * art * x(grid.index(i_r + 1, i_theta)); /* Top Right */
+                x(grid.index(i_r, i_theta_M1)) -= -0.25 * art * x(grid.index(i_r + 1, i_theta)); /* Top Right */
                 /* Fill result(i,j+1) */
-                x(grid.index(i_r, i_theta + 1)) -= +0.25 * art * x(grid.index(i_r + 1, i_theta)); /* Bottom Right */
+                x(grid.index(i_r, i_theta_P1)) -= +0.25 * art * x(grid.index(i_r + 1, i_theta)); /* Bottom Right */
 
                 /* -------------------------- */
                 /* Node on the outer boundary */
@@ -126,8 +133,8 @@ void DirectSolverGive<LevelCacheType>::applySymmetryShiftOuterBoundary(Vector<do
 
                 /* Fill result(i-1,j) */
                 x(grid.index(i_r - 1, i_theta)) -= -coeff1 * arr * x(grid.index(i_r, i_theta)) /* Right */
-                                                   - 0.25 * art * x(grid.index(i_r, i_theta + 1)) /* Top Right */
-                                                   + 0.25 * art * x(grid.index(i_r, i_theta - 1)); /* Bottom Right */
+                                                   - 0.25 * art * x(grid.index(i_r, i_theta_P1)) /* Top Right */
+                                                   + 0.25 * art * x(grid.index(i_r, i_theta_M1)); /* Bottom Right */
             }
         });
 }
diff --git a/include/GMGPolar/utils.h b/include/GMGPolar/utils.h
index df76f297..3823f80b 100644
--- a/include/GMGPolar/utils.h
+++ b/include/GMGPolar/utils.h
@@ -231,8 +231,9 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::writeToVTK(const std:
     file << "<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n";
     for (int i_r = 0; i_r < grid.nr() - 1; i_r++) {
         for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) {
+            const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
             file << grid.index(i_r, i_theta) << " " << grid.index(i_r + 1, i_theta) << " "
-                 << grid.index(i_r + 1, i_theta + 1) << " " << grid.index(i_r, i_theta + 1) << "\n";
+                 << grid.index(i_r + 1, i_theta_P1) << " " << grid.index(i_r, i_theta_P1) << "\n";
         }
     }
     file << "</DataArray>\n";
@@ -296,8 +297,9 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::writeToVTK(
     file << "<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n";
     for (int i_r = 0; i_r < grid.nr() - 1; i_r++) {
         for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) {
+            const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
             file << grid.index(i_r, i_theta) << " " << grid.index(i_r + 1, i_theta) << " "
-                 << grid.index(i_r + 1, i_theta + 1) << " " << grid.index(i_r, i_theta + 1) << "\n";
+                 << grid.index(i_r + 1, i_theta_P1) << " " << grid.index(i_r, i_theta_P1) << "\n";
         }
     }
 
diff --git a/include/PolarGrid/polargrid.h b/include/PolarGrid/polargrid.h
index 6f41feef..f6f0f4dc 100644
--- a/include/PolarGrid/polargrid.h
+++ b/include/PolarGrid/polargrid.h
@@ -45,7 +45,7 @@ class PolarGrid
 
     // Optimized, inlined indexing.
     KOKKOS_INLINE_FUNCTION int wrapThetaIndex(const int unwrapped_theta_index) const;
-    KOKKOS_INLINE_FUNCTION int index(const int r_index, const int unwrapped_theta_index) const;
+    KOKKOS_INLINE_FUNCTION int index(const int r_index, const int theta_index) const;
     KOKKOS_INLINE_FUNCTION void multiIndex(const int node_index, int& r_index, int& theta_index) const;
 
     // Grid Parameters
diff --git a/include/PolarGrid/polargrid.inl b/include/PolarGrid/polargrid.inl
index ada5fec0..309c643d 100644
--- a/include/PolarGrid/polargrid.inl
+++ b/include/PolarGrid/polargrid.inl
@@ -119,15 +119,16 @@ KOKKOS_INLINE_FUNCTION int PolarGrid::wrapThetaIndex(const int unwrapped_theta_i
     return theta_index;
 }
 
-KOKKOS_INLINE_FUNCTION int PolarGrid::index(const int r_index, const int unwrapped_theta_index) const
+KOKKOS_INLINE_FUNCTION int PolarGrid::index(const int r_index, const int theta_index) const
 {
-    // unwrapped_theta_index may be negative or larger than ntheta() to allow for periodicity.
     assert(0 <= r_index && r_index < nr());
-    int theta_index = wrapThetaIndex(unwrapped_theta_index);
+    assert(0 <= theta_index && theta_index < ntheta());
+
     int global_index =
         r_index < numberSmootherCircles()
             ? theta_index + ntheta() * r_index
             : numberCircularSmootherNodes() + r_index - numberSmootherCircles() + lengthRadialSmoother() * theta_index;
+
     assert(0 <= global_index && global_index < numberOfNodes());
     return global_index;
 }
diff --git a/src/Interpolation/extrapolated_prolongation.cpp b/src/Interpolation/extrapolated_prolongation.cpp
index c855bae6..64d01f80 100644
--- a/src/Interpolation/extrapolated_prolongation.cpp
+++ b/src/Interpolation/extrapolated_prolongation.cpp
@@ -56,11 +56,13 @@ static KOKKOS_INLINE_FUNCTION void fineNodeExtrapolatedProlongation(const int i_
     const int i_r_coarse     = i_r / 2;
     const int i_theta_coarse = i_theta / 2;
 
+    const int i_theta_coarse_P1 = coarse_grid.wrapThetaIndex(i_theta_coarse + 1);
+
     if (i_r & 1) {
         if (i_theta & 1) { /* (odd, odd) -> node in center of coarse cell */
             const double value =
                 0.5 * (coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse)] + /* Bottom right */
-                       coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] /* Top left */
+                       coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P1)] /* Top left */
                       );
             fine_result[fine_grid.index(i_r, i_theta)] = value;
         }
@@ -74,7 +76,7 @@ static KOKKOS_INLINE_FUNCTION void fineNodeExtrapolatedProlongation(const int i_
     else {
         if (i_theta & 1) { /* (even, odd) -> between coarse nodes in angular direction */
             const double value = 0.5 * (coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* Bottom */
-                                        coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] /* Top */
+                                        coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P1)] /* Top */
                                        );
             fine_result[fine_grid.index(i_r, i_theta)] = value;
         }
diff --git a/src/Interpolation/extrapolated_restriction.cpp b/src/Interpolation/extrapolated_restriction.cpp
index 5ed30ba6..6a497093 100644
--- a/src/Interpolation/extrapolated_restriction.cpp
+++ b/src/Interpolation/extrapolated_restriction.cpp
@@ -31,20 +31,23 @@ static KOKKOS_INLINE_FUNCTION void coarseNodeExtrapolatedRestriction(const int i
     const int i_r     = i_r_coarse * 2;
     const int i_theta = i_theta_coarse * 2;
 
+    const int i_theta_M1 = fine_grid.wrapThetaIndex(i_theta - 1);
+    const int i_theta_P1 = fine_grid.wrapThetaIndex(i_theta + 1);
+
     /* Center + Angular contributions (always present) */
-    double value = fine_values[fine_grid.index(i_r, i_theta)] + 0.5 * fine_values[fine_grid.index(i_r, i_theta - 1)] +
-                   0.5 * fine_values[fine_grid.index(i_r, i_theta + 1)];
+    double value = fine_values[fine_grid.index(i_r, i_theta)] + 0.5 * fine_values[fine_grid.index(i_r, i_theta_M1)] +
+                   0.5 * fine_values[fine_grid.index(i_r, i_theta_P1)];
 
     /* Left contributions (if not at inner boundary) */
     if (i_r_coarse > 0) {
         value += 0.5 * fine_values[fine_grid.index(i_r - 1, i_theta)] +
-                 0.5 * fine_values[fine_grid.index(i_r - 1, i_theta + 1)]; /* Top-Left diagonal */
+                 0.5 * fine_values[fine_grid.index(i_r - 1, i_theta_P1)]; /* Top-Left diagonal */
     }
 
     /* Right contributions (if not at outer boundary) */
     if (i_r_coarse < coarse_grid.nr() - 1) {
         value += 0.5 * fine_values[fine_grid.index(i_r + 1, i_theta)] +
-                 0.5 * fine_values[fine_grid.index(i_r + 1, i_theta - 1)]; /* Bottom-Right diagonal */
+                 0.5 * fine_values[fine_grid.index(i_r + 1, i_theta_M1)]; /* Bottom-Right diagonal */
     }
 
     coarse_result[coarse_grid.index(i_r_coarse, i_theta_coarse)] = value;
diff --git a/src/Interpolation/fmg_interpolation.cpp b/src/Interpolation/fmg_interpolation.cpp
index f9b3d174..5b0df52f 100644
--- a/src/Interpolation/fmg_interpolation.cpp
+++ b/src/Interpolation/fmg_interpolation.cpp
@@ -65,6 +65,10 @@ static KOKKOS_INLINE_FUNCTION void fineNodeFMGInterpolation(const int i_r, const
     const int i_r_coarse     = i_r / 2;
     const int i_theta_coarse = i_theta / 2;
 
+    const int i_theta_coarse_M1 = coarse_grid.wrapThetaIndex(i_theta_coarse - 1);
+    const int i_theta_coarse_P1 = coarse_grid.wrapThetaIndex(i_theta_coarse + 1);
+    const int i_theta_coarse_P2 = coarse_grid.wrapThetaIndex(i_theta_coarse + 2);
+
     /* Case 1: On the boundary */
     if (i_r == 0 || i_r == fine_grid.nr() - 1) {
         if (i_theta & 1) {
@@ -79,10 +83,10 @@ static KOKKOS_INLINE_FUNCTION void fineNodeFMGInterpolation(const int i_r, const
             const double w_theta3 = -(k0 + k1) / (k0 + k1 + k2 + k3) * k1 / (k1 + k2 + k3) * k2 / k3;
 
             fine_result[fine_grid.index(i_r, i_theta)] =
-                (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse - 1)] + /* (0, -3) */
+                (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_M1)] + /* (0, -3) */
                  w_theta1 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* (0, -1) */
-                 w_theta2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] + /* (0, +1) */
-                 w_theta3 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 2)] /* (0, +3) */
+                 w_theta2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P1)] + /* (0, +1) */
+                 w_theta3 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P2)] /* (0, +3) */
                 );
         }
         else {
@@ -103,20 +107,21 @@ static KOKKOS_INLINE_FUNCTION void fineNodeFMGInterpolation(const int i_r, const
             const double w_theta3 = -(k0 + k1) / (k0 + k1 + k2 + k3) * k1 / (k1 + k2 + k3) * k2 / k3;
 
             const double left_value =
-                (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse - 1)] + /* (-1, -3) */
+                (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_M1)] + /* (-1, -3) */
                  w_theta1 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* (-1, -1) */
-                 w_theta2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] + /* (-1, +1) */
-                 w_theta3 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 2)] /* (-1, +3) */
+                 w_theta2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P1)] + /* (-1, +1) */
+                 w_theta3 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P2)] /* (-1, +3) */
                 );
             const double right_value =
-                (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse - 1)] + /* (+1, -3) */
+                (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse_M1)] + /* (+1, -3) */
                  w_theta1 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse)] + /* (+1, -1) */
-                 w_theta2 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse + 1)] + /* (+1, +1) */
-                 w_theta3 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse + 2)] /* (+1, +3) */
+                 w_theta2 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse_P1)] + /* (+1, +1) */
+                 w_theta3 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse_P2)] /* (+1, +3) */
                 );
 
-            const double h1                            = fine_grid.radialSpacing(i_r - 1);
-            const double h2                            = fine_grid.radialSpacing(i_r);
+            const double h1 = fine_grid.radialSpacing(i_r - 1);
+            const double h2 = fine_grid.radialSpacing(i_r);
+
             fine_result[fine_grid.index(i_r, i_theta)] = (h1 * left_value + h2 * right_value) / (h1 + h2);
         }
         else {
@@ -143,28 +148,28 @@ static KOKKOS_INLINE_FUNCTION void fineNodeFMGInterpolation(const int i_r, const
                 const double w_theta3 = -(k0 + k1) / (k0 + k1 + k2 + k3) * k1 / (k1 + k2 + k3) * k2 / k3;
 
                 const double outer_left_value =
-                    (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse - 1, i_theta_coarse - 1)] + /* (-3, -3) */
+                    (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse - 1, i_theta_coarse_M1)] + /* (-3, -3) */
                      w_theta1 * coarse_values[coarse_grid.index(i_r_coarse - 1, i_theta_coarse)] + /* (-3, -1) */
-                     w_theta2 * coarse_values[coarse_grid.index(i_r_coarse - 1, i_theta_coarse + 1)] + /* (-3, +1) */
-                     w_theta3 * coarse_values[coarse_grid.index(i_r_coarse - 1, i_theta_coarse + 2)] /* (-3, +3) */
+                     w_theta2 * coarse_values[coarse_grid.index(i_r_coarse - 1, i_theta_coarse_P1)] + /* (-3, +1) */
+                     w_theta3 * coarse_values[coarse_grid.index(i_r_coarse - 1, i_theta_coarse_P2)] /* (-3, +3) */
                     );
                 const double inner_left_value =
-                    (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse - 1)] + /* (-1, -3) */
+                    (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_M1)] + /* (-1, -3) */
                      w_theta1 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* (-1, -1) */
-                     w_theta2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] + /* (-1, +1) */
-                     w_theta3 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 2)] /* (-1, +3) */
+                     w_theta2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P1)] + /* (-1, +1) */
+                     w_theta3 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P2)] /* (-1, +3) */
                     );
                 const double inner_right_value =
-                    (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse - 1)] + /* (+1, -3) */
+                    (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse_M1)] + /* (+1, -3) */
                      w_theta1 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse)] + /* (+1, -1) */
-                     w_theta2 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse + 1)] + /* (+1, +1) */
-                     w_theta3 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse + 2)] /* (+1, +3) */
+                     w_theta2 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse_P1)] + /* (+1, +1) */
+                     w_theta3 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse_P2)] /* (+1, +3) */
                     );
                 const double outer_right_value =
-                    (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse + 2, i_theta_coarse - 1)] + /* (+3, -3) */
+                    (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse + 2, i_theta_coarse_M1)] + /* (+3, -3) */
                      w_theta1 * coarse_values[coarse_grid.index(i_r_coarse + 2, i_theta_coarse)] + /* (+3, -1) */
-                     w_theta2 * coarse_values[coarse_grid.index(i_r_coarse + 2, i_theta_coarse + 1)] + /* (+3, +1) */
-                     w_theta3 * coarse_values[coarse_grid.index(i_r_coarse + 2, i_theta_coarse + 2)] /* (+3, +3) */
+                     w_theta2 * coarse_values[coarse_grid.index(i_r_coarse + 2, i_theta_coarse_P1)] + /* (+3, +1) */
+                     w_theta3 * coarse_values[coarse_grid.index(i_r_coarse + 2, i_theta_coarse_P2)] /* (+3, +3) */
                     );
 
                 const double h0 = coarse_grid.radialSpacing(i_r_coarse - 1);
@@ -212,10 +217,10 @@ static KOKKOS_INLINE_FUNCTION void fineNodeFMGInterpolation(const int i_r, const
                 const double w_theta3 = -(k0 + k1) / (k0 + k1 + k2 + k3) * k1 / (k1 + k2 + k3) * k2 / k3;
 
                 fine_result[fine_grid.index(i_r, i_theta)] =
-                    (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse - 1)] + /* (0, -3) */
+                    (w_theta0 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_M1)] + /* (0, -3) */
                      w_theta1 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* (0, -1) */
-                     w_theta2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] + /* (0, +1) */
-                     w_theta3 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 2)] /* (0, +3) */
+                     w_theta2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P1)] + /* (0, +1) */
+                     w_theta3 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P2)] /* (0, +3) */
                     );
             }
             else {
diff --git a/src/Interpolation/prolongation.cpp b/src/Interpolation/prolongation.cpp
index fc7fd203..441d7d1f 100644
--- a/src/Interpolation/prolongation.cpp
+++ b/src/Interpolation/prolongation.cpp
@@ -59,6 +59,8 @@ static KOKKOS_INLINE_FUNCTION void fineNodeProlongation(const int i_r, const int
     const int i_r_coarse     = i_r / 2;
     const int i_theta_coarse = i_theta / 2;
 
+    const int i_theta_coarse_P1 = coarse_grid.wrapThetaIndex(i_theta_coarse + 1);
+
     if (i_r & 1) {
         if (i_theta & 1) { /* (odd, odd) -> fine node in center of coarse cell */
             const double h1 = fine_grid.radialSpacing(i_r - 1);
@@ -69,8 +71,8 @@ static KOKKOS_INLINE_FUNCTION void fineNodeProlongation(const int i_r, const int
             const double value =
                 (h1 * k1 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* Bottom left */
                  h2 * k1 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse)] + /* Bottom right */
-                 h1 * k2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] + /* Top left */
-                 h2 * k2 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse + 1)] /* Top right */
+                 h1 * k2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P1)] + /* Top left */
+                 h2 * k2 * coarse_values[coarse_grid.index(i_r_coarse + 1, i_theta_coarse_P1)] /* Top right */
                  ) /
                 ((h1 + h2) * (k1 + k2));
 
@@ -94,7 +96,7 @@ static KOKKOS_INLINE_FUNCTION void fineNodeProlongation(const int i_r, const int
             const double k2 = fine_grid.angularSpacing(i_theta);
 
             const double value = (k1 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse)] + /* Bottom */
-                                  k2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse + 1)] /* Top */
+                                  k2 * coarse_values[coarse_grid.index(i_r_coarse, i_theta_coarse_P1)] /* Top */
                                   ) /
                                  (k1 + k2);
 
diff --git a/tests/GMGPolar/pcg_tests.cpp b/tests/GMGPolar/pcg_tests.cpp
index 2c3a0ccf..3b175afb 100644
--- a/tests/GMGPolar/pcg_tests.cpp
+++ b/tests/GMGPolar/pcg_tests.cpp
@@ -196,7 +196,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.5>, // refinementRadius
         std::integral_constant<int, 1>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
@@ -238,7 +238,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.5>, // refinementRadius
         std::integral_constant<int, 1>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -280,7 +280,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.5>, // refinementRadius
         std::integral_constant<int, 1>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_TAKE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -320,7 +320,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.6>, // refinementRadius
         std::integral_constant<int, 3>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_TAKE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -360,7 +360,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66>, // refinementRadius
         std::integral_constant<int, 0>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -400,7 +400,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66>, // refinementRadius
         std::integral_constant<int, 1>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
@@ -440,7 +440,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66>, // refinementRadius
         std::integral_constant<int, 1>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
@@ -479,7 +479,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66>, // refinementRadius
         std::integral_constant<int, 1>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -518,7 +518,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.6>, // refinementRadius
         std::integral_constant<int, 2>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_TAKE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -557,7 +557,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66>, // refinementRadius
         std::integral_constant<int, 2>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -596,7 +596,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66 * 1.2>, // refinementRadius
         std::integral_constant<int, 0>, // anisotropicFactor
         std::integral_constant<int, 1>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
diff --git a/tests/GMGPolar/solve_tests.cpp b/tests/GMGPolar/solve_tests.cpp
index 074bb340..ebe83db8 100644
--- a/tests/GMGPolar/solve_tests.cpp
+++ b/tests/GMGPolar/solve_tests.cpp
@@ -172,7 +172,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.5>, // refinementRadius
         std::integral_constant<int, 1>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
@@ -208,7 +208,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.5>, // refinementRadius
         std::integral_constant<int, 1>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -244,7 +244,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.5>, // refinementRadius
         std::integral_constant<int, 1>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_TAKE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -278,7 +278,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.6>, // refinementRadius
         std::integral_constant<int, 3>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_TAKE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -312,7 +312,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66>, // refinementRadius
         std::integral_constant<int, 0>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -346,7 +346,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66>, // refinementRadius
         std::integral_constant<int, 1>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
@@ -380,7 +380,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66>, // refinementRadius
         std::integral_constant<int, 1>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
@@ -413,7 +413,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66>, // refinementRadius
         std::integral_constant<int, 1>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -446,7 +446,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.6>, // refinementRadius
         std::integral_constant<int, 2>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_TAKE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -479,7 +479,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66>, // refinementRadius
         std::integral_constant<int, 2>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -512,7 +512,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66 * 1.2>, // refinementRadius
         std::integral_constant<int, 0>, // anisotropicFactor
         std::integral_constant<int, 1>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -545,7 +545,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66 * 1.2>, // refinementRadius
         std::integral_constant<int, 0>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -578,7 +578,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66 * 1.2>, // refinementRadius
         std::integral_constant<int, 0>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -611,7 +611,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 0.66 * 1.2>, // refinementRadius
         std::integral_constant<int, 0>, // anisotropicFactor
         std::integral_constant<int, 0>, // divideBy2
-        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
         std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
@@ -654,7 +654,7 @@ void run_gmgpolar()
 
     GMGPolar solver(grid, domain, profile_coefficients);
 
-    bool paraview          = false;
+    bool paraview          = true;
     int preSmoothingSteps  = 1;
     int postSmoothingSteps = 1;
 
diff --git a/tests/Interpolation/extrapolated_prolongation.cpp b/tests/Interpolation/extrapolated_prolongation.cpp
index 160178d4..4ee51d38 100644
--- a/tests/Interpolation/extrapolated_prolongation.cpp
+++ b/tests/Interpolation/extrapolated_prolongation.cpp
@@ -18,6 +18,8 @@ static double expected_extrapolated_value(const PolarGrid& coarse, const PolarGr
     bool r_even = (i_r % 2 == 0);
     bool t_even = (i_theta % 2 == 0);
 
+    int i_theta_coarse_P1 = coarse.wrapThetaIndex(i_theta_coarse + 1);
+
     double result = 0;
     Kokkos::parallel_reduce(
         "extrap_prolongation_test", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, 1),
@@ -36,12 +38,12 @@ static double expected_extrapolated_value(const PolarGrid& coarse, const PolarGr
             else if (r_even && !t_even) {
                 // Angular midpoint - arithmetic mean of bottom and top
                 local_result = 0.5 * (coarse_vals[coarse.index(i_r_coarse, i_theta_coarse)] +
-                                      coarse_vals[coarse.index(i_r_coarse, i_theta_coarse + 1)]);
+                                      coarse_vals[coarse.index(i_r_coarse, i_theta_coarse_P1)]);
             }
             else {
                 // Center of coarse cell - arithmetic mean of diagonal nodes (bottom-right + top-left)
                 local_result = 0.5 * (coarse_vals[coarse.index(i_r_coarse + 1, i_theta_coarse)] +
-                                      coarse_vals[coarse.index(i_r_coarse, i_theta_coarse + 1)]);
+                                      coarse_vals[coarse.index(i_r_coarse, i_theta_coarse_P1)]);
             }
         },
         result);
diff --git a/tests/Interpolation/prolongation.cpp b/tests/Interpolation/prolongation.cpp
index ccff5740..44335c65 100644
--- a/tests/Interpolation/prolongation.cpp
+++ b/tests/Interpolation/prolongation.cpp
@@ -17,6 +17,8 @@ static double expected_value(const PolarGrid& coarse, const PolarGrid& fine, Con
     bool r_even = (i_r % 2 == 0);
     bool t_even = (i_theta % 2 == 0);
 
+    int i_theta_coarse_P1 = coarse.wrapThetaIndex(i_theta_coarse + 1);
+
     double result = 0;
     Kokkos::parallel_reduce(
         "prolongation_test", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, 1),
@@ -41,8 +43,8 @@ static double expected_value(const PolarGrid& coarse, const PolarGrid& fine, Con
                 double k1 = fine.angularSpacing(i_theta - 1);
                 double k2 = fine.angularSpacing(i_theta);
 
-                local_result = (k1 * coarse_vals[coarse.index(i_r_coarse, t_even ? i_theta_coarse : i_theta_coarse)] +
-                                k2 * coarse_vals[coarse.index(i_r_coarse, i_theta_coarse + 1)]) /
+                local_result = (k1 * coarse_vals[coarse.index(i_r_coarse, i_theta_coarse)] +
+                                k2 * coarse_vals[coarse.index(i_r_coarse, i_theta_coarse_P1)]) /
                                (k1 + k2);
             }
             else {
@@ -54,8 +56,8 @@ static double expected_value(const PolarGrid& coarse, const PolarGrid& fine, Con
 
                 local_result = (h1 * k1 * coarse_vals[coarse.index(i_r_coarse, i_theta_coarse)] +
                                 h2 * k1 * coarse_vals[coarse.index(i_r_coarse + 1, i_theta_coarse)] +
-                                h1 * k2 * coarse_vals[coarse.index(i_r_coarse, i_theta_coarse + 1)] +
-                                h2 * k2 * coarse_vals[coarse.index(i_r_coarse + 1, i_theta_coarse + 1)]) /
+                                h1 * k2 * coarse_vals[coarse.index(i_r_coarse, i_theta_coarse_P1)] +
+                                h2 * k2 * coarse_vals[coarse.index(i_r_coarse + 1, i_theta_coarse_P1)]) /
                                ((h1 + h2) * (k1 + k2));
             }
         },