diff --git a/include/Definitions/geometry_helper.h b/include/Definitions/geometry_helper.h
index 1245e598..67aa739b 100644
--- a/include/Definitions/geometry_helper.h
+++ b/include/Definitions/geometry_helper.h
@@ -33,4 +33,7 @@ KOKKOS_INLINE_FUNCTION void compute_jacobian_elements(const DomainGeometry& doma
     /*  1.0 / det(DF) *  */
     /*  [Jtt, -Jrt]      */
     /*  [-Jtr, Jrr]      */
+    /* The stencil coefficient art is scaled by 0.25 to account for the */
+    /* 9-point stencil used in the finite difference approximation */
+    art *= 0.25;
 }
diff --git a/include/DirectSolver/DirectSolverGive/applySymmetryShift.inl b/include/DirectSolver/DirectSolverGive/applySymmetryShift.inl
index 56f48a9e..13aa872d 100644
--- a/include/DirectSolver/DirectSolverGive/applySymmetryShift.inl
+++ b/include/DirectSolver/DirectSolverGive/applySymmetryShift.inl
@@ -45,8 +45,8 @@ void DirectSolverGive<LevelCacheType>::applySymmetryShiftInnerBoundary(Vector<do
 
                 /* 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_P1)) /* Top Left */
-                                                   - 0.25 * art * x(grid.index(i_r, i_theta_M1)); /* Bottom Left */
+                                                   + art * x(grid.index(i_r, i_theta_P1)) /* Top Left */
+                                                   - art * x(grid.index(i_r, i_theta_M1)); /* Bottom Left */
 
                 /* --------------------------- */
                 /* Node next to inner boundary */
@@ -66,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_M1)) -= +0.25 * art * x(grid.index(i_r - 1, i_theta)); /* Top Left */
+                x(grid.index(i_r, i_theta_M1)) -= +art * x(grid.index(i_r - 1, i_theta)); /* Top Left */
                 /* Fill x(i,j+1) */
-                x(grid.index(i_r, i_theta_P1)) -= -0.25 * art * x(grid.index(i_r - 1, i_theta)); /* Bottom Left */
+                x(grid.index(i_r, i_theta_P1)) -= -art * x(grid.index(i_r - 1, i_theta)); /* Bottom Left */
             }
         });
 }
@@ -112,9 +112,9 @@ void DirectSolverGive<LevelCacheType>::applySymmetryShiftOuterBoundary(Vector<do
                 /* 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_M1)) -= -0.25 * art * x(grid.index(i_r + 1, i_theta)); /* Top Right */
+                x(grid.index(i_r, i_theta_M1)) -= -art * x(grid.index(i_r + 1, i_theta)); /* Top Right */
                 /* Fill result(i,j+1) */
-                x(grid.index(i_r, i_theta_P1)) -= +0.25 * art * x(grid.index(i_r + 1, i_theta)); /* Bottom Right */
+                x(grid.index(i_r, i_theta_P1)) -= +art * x(grid.index(i_r + 1, i_theta)); /* Bottom Right */
 
                 /* -------------------------- */
                 /* Node on the outer boundary */
@@ -133,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_P1)) /* Top Right */
-                                                   + 0.25 * art * x(grid.index(i_r, i_theta_M1)); /* Bottom Right */
+                                                   - art * x(grid.index(i_r, i_theta_P1)) /* Top Right */
+                                                   + art * x(grid.index(i_r, i_theta_M1)); /* Bottom Right */
             }
         });
 }
diff --git a/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl b/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl
index a1c02bea..18969920 100644
--- a/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl
+++ b/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl
@@ -126,12 +126,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
         offset = LeftStencil[StencilPosition::TopRight];
         column = top_index;
-        value  = -0.25 * art; /* Top Right */
+        value  = -art; /* Top Right */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         offset = LeftStencil[StencilPosition::BottomRight];
         column = bottom_index;
-        value  = +0.25 * art; /* Bottom Right */
+        value  = +art; /* Bottom Right */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         /* Fill matrix row of (i+1,j) */
@@ -152,12 +152,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
         offset = RightStencil[StencilPosition::TopLeft];
         column = top_index;
-        value  = +0.25 * art; /* Top Left */
+        value  = +art; /* Top Left */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         offset = RightStencil[StencilPosition::BottomLeft];
         column = bottom_index;
-        value  = -0.25 * art; /* Bottom Left */
+        value  = -art; /* Bottom Left */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         /* Fill matrix row of (i,j-1) */
@@ -178,12 +178,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
         offset = BottomStencil[StencilPosition::TopRight];
         column = right_index;
-        value  = -0.25 * art; /* Top Right */
+        value  = -art; /* Top Right */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         offset = BottomStencil[StencilPosition::TopLeft];
         column = left_index;
-        value  = +0.25 * art; /* Top Left */
+        value  = +art; /* Top Left */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         /* Fill matrix row of (i,j+1) */
@@ -204,12 +204,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
         offset = TopStencil[StencilPosition::BottomRight];
         column = right_index;
-        value  = +0.25 * art; /* Bottom Right */
+        value  = +art; /* Bottom Right */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         offset = TopStencil[StencilPosition::BottomLeft];
         column = left_index;
-        value  = -0.25 * art; /* Bottom Left */
+        value  = -art; /* Bottom Left */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
     }
     /* -------------------------- */
@@ -369,12 +369,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
             offset = RightStencil[StencilPosition::TopLeft];
             column = top_index;
-            value  = +0.25 * art; /* Top Left */
+            value  = +art; /* Top Left */
             updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
             offset = RightStencil[StencilPosition::BottomLeft];
             column = bottom_index;
-            value  = -0.25 * art; /* Bottom Left */
+            value  = -art; /* Bottom Left */
             updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
             /* Fill matrix row of (i,j-1) */
@@ -395,7 +395,7 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
             offset = BottomStencil[StencilPosition::TopRight];
             column = right_index;
-            value  = -0.25 * art; /* Top Right */
+            value  = -art; /* Top Right */
             updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
             /* TopLeft REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
@@ -418,7 +418,7 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
             offset = TopStencil[StencilPosition::BottomRight];
             column = right_index;
-            value  = +0.25 * art; /* Bottom Right */
+            value  = +art; /* Bottom Right */
             updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
             /* BottomLeft REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
@@ -512,12 +512,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
             offset = LeftStencil[StencilPosition::TopRight];
             column = top_index;
-            value  = -0.25 * art; /* Top Right */
+            value  = -art; /* Top Right */
             updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
             offset = LeftStencil[StencilPosition::BottomRight];
             column = bottom_index;
-            value  = +0.25 * art; /* Bottom Right */
+            value  = +art; /* Bottom Right */
             updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
         }
 
@@ -539,12 +539,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
         offset = RightStencil[StencilPosition::TopLeft];
         column = top_index;
-        value  = +0.25 * art; /* Top Left */
+        value  = +art; /* Top Left */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         offset = RightStencil[StencilPosition::BottomLeft];
         column = bottom_index;
-        value  = -0.25 * art; /* Bottom Left */
+        value  = -art; /* Bottom Left */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         /* Fill matrix row of (i,j-1) */
@@ -565,14 +565,14 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
         offset = BottomStencil[StencilPosition::TopRight];
         column = right_index;
-        value  = -0.25 * art; /* Top Right */
+        value  = -art; /* Top Right */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         /* REMOVED: Moved to the right hand side to make the matrix symmetric */
         if (!DirBC_Interior) {
             offset = BottomStencil[StencilPosition::TopLeft];
             column = left_index;
-            value  = +0.25 * art; /* Top Left */
+            value  = +art; /* Top Left */
             updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
         }
 
@@ -594,14 +594,14 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
         offset = TopStencil[StencilPosition::BottomRight];
         column = right_index;
-        value  = +0.25 * art; /* Bottom Right */
+        value  = +art; /* Bottom Right */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         /* REMOVED: Moved to the right hand side to make the matrix symmetric */
         if (!DirBC_Interior) {
             offset = TopStencil[StencilPosition::BottomLeft];
             column = left_index;
-            value  = -0.25 * art; /* Bottom Left */
+            value  = -art; /* Bottom Left */
             updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
         }
     }
@@ -684,12 +684,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
         offset = LeftStencil[StencilPosition::TopRight];
         column = top_index;
-        value  = -0.25 * art; /* Top Right */
+        value  = -art; /* Top Right */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         offset = LeftStencil[StencilPosition::BottomRight];
         column = bottom_index;
-        value  = 0.25 * art; /* Bottom Right */
+        value  = art; /* Bottom Right */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         /* Fill matrix row of (i+1,j) */
@@ -715,7 +715,7 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
         offset = BottomStencil[StencilPosition::TopLeft];
         column = left_index;
-        value  = 0.25 * art; /* Top Left */
+        value  = art; /* Top Left */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
 
         /* Fill matrix row of (i,j+1) */
@@ -738,7 +738,7 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri
 
         offset = TopStencil[StencilPosition::BottomLeft];
         column = left_index;
-        value  = -0.25 * art; /* Bottom Left */
+        value  = -art; /* Bottom Left */
         updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
     }
     /* ------------------------------------ */
diff --git a/include/DirectSolver/DirectSolverTake/applySymmetryShift.inl b/include/DirectSolver/DirectSolverTake/applySymmetryShift.inl
index 775febff..84e76a6b 100644
--- a/include/DirectSolver/DirectSolverTake/applySymmetryShift.inl
+++ b/include/DirectSolver/DirectSolverTake/applySymmetryShift.inl
@@ -43,8 +43,8 @@ void DirectSolverTake<LevelCacheType>::applySymmetryShiftInnerBoundary(Vector<do
             const int top    = grid.index(i_r, i_theta_P1);
 
             x[center] -= (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */
-                          - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
-                          + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */
+                          - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
+                          + (art[left] + art[top]) * x[top_left] /* Top Left */
             );
         });
 }
@@ -86,8 +86,8 @@ void DirectSolverTake<LevelCacheType>::applySymmetryShiftOuterBoundary(Vector<do
             const int top_right    = grid.index(i_r + 1, i_theta_P1);
 
             x[center] -= (-coeff2 * (arr[center] + arr[right]) * x[right] /* Right */
-                          + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
-                          - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */
+                          + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
+                          - (art[right] + art[top]) * x[top_right] /* Top Right */
             );
         });
 }
diff --git a/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl b/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl
index 3f04c610..058c3ef1 100644
--- a/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl
+++ b/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl
@@ -78,10 +78,10 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri
              - top_value /* Center: (Top) */
             );
 
-        const double bottom_left_value  = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */
-        const double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */
-        const double top_left_value     = +0.25 * (art(left_index) + art(top_index)); /* Top Left */
-        const double top_right_value    = -0.25 * (art(right_index) + art(top_index)); /* Top Right */
+        const double bottom_left_value  = -(art(left_index) + art(bottom_index)); /* Bottom Left */
+        const double bottom_right_value = +(art(right_index) + art(bottom_index)); /* Bottom Right */
+        const double top_left_value     = +(art(left_index) + art(top_index)); /* Top Left */
+        const double top_right_value    = -(art(right_index) + art(top_index)); /* Top Right */
 
         /* Fill matrix row of (i,j) */
         row = center_index;
@@ -204,8 +204,8 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri
                  - top_value /* Center: (Top) */
                 );
 
-            const double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */
-            const double top_right_value    = -0.25 * (art(right_index) + art(top_index)); /* Top Right */
+            const double bottom_right_value = +(art(right_index) + art(bottom_index)); /* Bottom Right */
+            const double top_right_value    = -(art(right_index) + art(top_index)); /* Top Right */
 
             /* Fill matrix row of (i,j) */
             row = center_index;
@@ -296,10 +296,10 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri
              - top_value /* Center: (Top) */
             );
 
-        const double bottom_left_value  = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */
-        const double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */
-        const double top_left_value     = +0.25 * (art(left_index) + art(top_index)); /* Top Left */
-        const double top_right_value    = -0.25 * (art(right_index) + art(top_index)); /* Top Right */
+        const double bottom_left_value  = -(art(left_index) + art(bottom_index)); /* Bottom Left */
+        const double bottom_right_value = +(art(right_index) + art(bottom_index)); /* Bottom Right */
+        const double top_left_value     = +(art(left_index) + art(top_index)); /* Top Left */
+        const double top_right_value    = -(art(right_index) + art(top_index)); /* Top Right */
 
         /* Fill matrix row of (i,j) */
         row = center_index;
@@ -402,8 +402,8 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri
              - top_value /* Center: (Top) */
             );
 
-        const double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */
-        const double top_left_value    = +0.25 * (art(left_index) + art(top_index)); /* Top Left */
+        const double bottom_left_value = -(art(left_index) + art(bottom_index)); /* Bottom Left */
+        const double top_left_value    = +(art(left_index) + art(top_index)); /* Top Left */
 
         /* Fill matrix row of (i,j) */
         row = center_index;
diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/applyAscOrtho.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/applyAscOrtho.inl
index ae389c24..5598f2e2 100644
--- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/applyAscOrtho.inl
+++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/applyAscOrtho.inl
@@ -62,11 +62,11 @@ nodeApplyAscOrthoCircleGiveInside(const int i_r, const int i_theta, const PolarG
             result[center] -= (-coeff1 * arr * x[left] /* Left */
                                - coeff2 * arr * x[right] /* Right */
             ); /* Fill result(i,j-1) */
-            result[bottom] -= (-0.25 * art * x[right] /* Top Right */
-                               + 0.25 * art * x[left]); /* Top Left */
+            result[bottom] -= (-art * x[right] /* Top Right */
+                               + art * x[left]); /* Top Left */
             /* Fill result(i,j+1) */
-            result[top] -= (+0.25 * art * x[right] /* Bottom Right */
-                            - 0.25 * art * x[left]); /* Bottom Left */
+            result[top] -= (+art * x[right] /* Bottom Right */
+                            - art * x[left]); /* Bottom Left */
         }
         else {
             if (i_theta & 1) {
@@ -94,12 +94,12 @@ nodeApplyAscOrthoCircleGiveInside(const int i_r, const int i_theta, const PolarG
 
                 /* Fill result(i,j-1) */
                 result[bottom] -= (-coeff3 * att * x[center] /* Top */
-                                   - 0.25 * art * x[right] /* Top Right */
-                                   + 0.25 * art * x[left]); /* Top Left */
+                                   - art * x[right] /* Top Right */
+                                   + art * x[left]); /* Top Left */
                 /* Fill result(i,j+1) */
                 result[top] -= (-coeff4 * att * x[center] /* Bottom */
-                                + 0.25 * art * x[right] /* Bottom Right */
-                                - 0.25 * art * x[left]); /* Bottom Left */
+                                + art * x[right] /* Bottom Right */
+                                - art * x[left]); /* Bottom Left */
             }
         }
     }
@@ -184,12 +184,12 @@ nodeApplyAscOrthoCircleGiveInside(const int i_r, const int i_theta, const PolarG
 
                 /* Fill result(i,j-1) */
                 result[bottom] -= (-coeff3 * att * x[center] /* Top */
-                                   - 0.25 * art * x[right]); /* Top Right */
-                /*  + 0.25 * art * x[grid.index(i_r-1,i_theta)]; // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
+                                   - art * x[right]); /* Top Right */
+                /*  + art * x[grid.index(i_r-1,i_theta)]; // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
                 /* Fill result(i,j+1) */
                 result[top] -= (-coeff4 * att * x[center] /* Bottom */
-                                + 0.25 * art * x[right]); /* Bottom Right */
-                /*  - 0.25 * art * x[grid.index(i_r-1,i_theta)]; // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
+                                + art * x[right]); /* Bottom Right */
+                /*  - art * x[grid.index(i_r-1,i_theta)]; // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
             }
         }
     }
@@ -346,8 +346,8 @@ nodeApplyAscOrthoCircleGiveOutside(const int i_r, const int i_theta, const Polar
 
         /* Fill result(i-1,j) */
         result[left] -= (-coeff1 * arr * x[center] /* Right */
-                         - 0.25 * art * x[top] /* Top Right */
-                         + 0.25 * art * x[bottom]); /* Bottom Right */
+                         - art * x[top] /* Top Right */
+                         + art * x[bottom]); /* Bottom Right */
     }
 
     if (give_right) {
@@ -357,8 +357,8 @@ nodeApplyAscOrthoCircleGiveOutside(const int i_r, const int i_theta, const Polar
 
         /* Fill result(i+1,j) */
         result[right] -= (-coeff2 * arr * x[center] /* Left */
-                          + 0.25 * art * x[top] /* Top Left */
-                          - 0.25 * art * x[bottom]); /* Bottom Left */
+                          + art * x[top] /* Top Left */
+                          - art * x[bottom]); /* Bottom Left */
     }
 }
 
@@ -426,11 +426,11 @@ nodeApplyAscOrthoRadialGiveInside(const int i_r, const int i_theta, const PolarG
             result[center] -= (-coeff3 * att * x[bottom] /* Bottom */
                                - coeff4 * att * x[top]); /* Top */
             /* Fill result(i-1,j) */
-            result[left] -= (-0.25 * art * x[top] /* Top Right */
-                             + 0.25 * art * x[bottom]); /* Bottom Right */
+            result[left] -= (-art * x[top] /* Top Right */
+                             + art * x[bottom]); /* Bottom Right */
             /* Fill result(i+1,j) */
-            result[right] -= (+0.25 * art * x[top] /* Top Left */
-                              - 0.25 * art * x[bottom]); /* Bottom Left */
+            result[right] -= (+art * x[top] /* Top Left */
+                              - art * x[bottom]); /* Bottom Left */
         }
         else {
             if (i_r & 1) {
@@ -461,12 +461,12 @@ nodeApplyAscOrthoRadialGiveInside(const int i_r, const int i_theta, const PolarG
 
                 /* Fill result(i-1,j) */
                 result[left] -= (-coeff1 * arr * x[center] /* Right */
-                                 - 0.25 * art * x[top] /* Top Right */
-                                 + 0.25 * art * x[bottom]); /* Bottom Right */
+                                 - art * x[top] /* Top Right */
+                                 + art * x[bottom]); /* Bottom Right */
                 /* Fill result(i+1,j) */
                 result[right] -= (-coeff2 * arr * x[center] /* Left */
-                                  + 0.25 * art * x[top] /* Top Left */
-                                  - 0.25 * art * x[bottom]); /* Bottom Left */
+                                  + art * x[top] /* Top Left */
+                                  - art * x[bottom]); /* Bottom Left */
             }
         }
     }
@@ -494,8 +494,8 @@ nodeApplyAscOrthoRadialGiveInside(const int i_r, const int i_theta, const PolarG
         if ((!(i_r & 1) || (i_theta & 1))) {
             /* Fill result(i+1,j) */
             result[right] -= (-coeff2 * arr * x[center] /* Left */
-                              + 0.25 * art * x[top] /* Top Left */
-                              - 0.25 * art * x[bottom]); /* Bottom Left */
+                              + art * x[top] /* Top Left */
+                              - art * x[bottom]); /* Bottom Left */
         }
     }
     else if (i_r == grid.numberSmootherCircles()) {
@@ -532,8 +532,8 @@ nodeApplyAscOrthoRadialGiveInside(const int i_r, const int i_theta, const PolarG
                                    - coeff4 * att * x[top] /* Top */
                 );
                 /* Fill result(i+1,j) */
-                result[right] -= (+0.25 * art * x[top] /* Top Left */
-                                  - 0.25 * art * x[bottom]); /* Bottom Left */
+                result[right] -= (+art * x[top] /* Top Left */
+                                  - art * x[bottom]); /* Bottom Left */
             }
             else {
                 /* i_theta % 2 == 1 and i_r % 2 == 0 */
@@ -549,8 +549,8 @@ nodeApplyAscOrthoRadialGiveInside(const int i_r, const int i_theta, const PolarG
                                    - coeff4 * att * x[top] /* Top */
                 );
                 /* Fill result(i+1,j) */
-                result[right] -= (+0.25 * art * x[top] /* Top Left */
-                                  - 0.25 * art * x[bottom]); /* Bottom Left */
+                result[right] -= (+art * x[top] /* Top Left */
+                                  - art * x[bottom]); /* Bottom Left */
             }
         }
         else {
@@ -579,8 +579,8 @@ nodeApplyAscOrthoRadialGiveInside(const int i_r, const int i_theta, const PolarG
 
                 /* Fill result(i+1,j) */
                 result[right] -= (-coeff2 * arr * x[center] /* Left */
-                                  + 0.25 * art * x[top] /* Top Left */
-                                  - 0.25 * art * x[bottom]); /* Bottom Left */
+                                  + art * x[top] /* Top Left */
+                                  - art * x[bottom]); /* Bottom Left */
             }
         }
     }
@@ -620,8 +620,8 @@ nodeApplyAscOrthoRadialGiveInside(const int i_r, const int i_theta, const PolarG
                                - coeff4 * att * x[top] /* Top */
             );
             /* Fill result(i-1,j) */
-            result[left] -= (-0.25 * art * x[top] /* Top Right */
-                             + 0.25 * art * x[bottom]); /* Bottom Right */
+            result[left] -= (-art * x[top] /* Top Right */
+                             + art * x[bottom]); /* Bottom Right */
 
             /* "Right" is part of the radial Asc smoother matrices, */
             /* but is shifted over to the rhs to make the radial Asc smoother matrices symmetric. */
@@ -673,8 +673,8 @@ nodeApplyAscOrthoRadialGiveInside(const int i_r, const int i_theta, const PolarG
             /* -----------|| */
 
             /* Fill result(i-1,j) */
-            result[left] -= (-0.25 * art * x[top] /* Top Right */
-                             + 0.25 * art * x[bottom] /* Bottom Right */
+            result[left] -= (-art * x[top] /* Top Right */
+                             + art * x[bottom] /* Bottom Right */
             );
             /* "Right" is part of the radial Asc smoother matrices, */
             /* but is shifted over to the rhs to make the radial Asc smoother matrices symmetric. */
@@ -694,8 +694,8 @@ nodeApplyAscOrthoRadialGiveInside(const int i_r, const int i_theta, const PolarG
 
             /* Fill result(i-1,j) */
             result[left] -= (-coeff1 * arr * x[center] /* Right */
-                             - 0.25 * art * x[top] /* Top Right */
-                             + 0.25 * art * x[bottom]); /* Bottom Right */
+                             - art * x[top] /* Top Right */
+                             + art * x[bottom]); /* Bottom Right */
         }
     }
 }
@@ -821,8 +821,8 @@ nodeApplyAscOrthoRadialGiveOutside(const int i_r, const int i_theta, const Polar
 
         /* Fill result(i,j-1) */
         result[bottom] -= (-coeff3 * att * x[center] /* Top */
-                           - 0.25 * art * x[right] /* Top Right */
-                           + 0.25 * art * x[left]); /* Top Left */
+                           - art * x[right] /* Top Right */
+                           + art * x[left]); /* Top Left */
     }
 
     if (give_top) {
@@ -831,8 +831,8 @@ nodeApplyAscOrthoRadialGiveOutside(const int i_r, const int i_theta, const Polar
         const int top       = grid.index(i_r, i_theta_P1);
         /* Fill result(i,j+1) */
         result[top] -= (-coeff4 * att * x[center] /* Bottom */
-                        + 0.25 * art * x[right] /* Bottom Right */
-                        - 0.25 * art * x[left]); /* Bottom Left */
+                        + art * x[right] /* Bottom Right */
+                        - art * x[left]); /* Bottom Left */
     }
 }
 
diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl
index ef66592f..da709fc4 100644
--- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl
+++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl
@@ -3,205 +3,86 @@
 namespace extrapolated_smoother_take
 {
 
-static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleTake(const int i_r, const int i_theta, const PolarGrid& grid,
-                                                               const bool DirBC_Interior, ConstVector<double>& x,
-                                                               ConstVector<double>& rhs, Vector<double>& result,
-                                                               ConstVector<double>& arr, ConstVector<double>& att,
-                                                               ConstVector<double>& art)
+static KOKKOS_INLINE_FUNCTION void
+nodeApplyAscOrthoCircleTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, ConstVector<double>& x,
+                                    ConstVector<double>& rhs, Vector<double>& result, ConstVector<double>& arr,
+                                    ConstVector<double>& att, ConstVector<double>& art)
 {
-    KOKKOS_ASSERT(i_r >= 0 && i_r <= grid.numberSmootherCircles());
-
-    /* -------------------- */
-    /* Node in the interior */
-    /* -------------------- */
-    if (i_r > 0 && i_r < grid.numberSmootherCircles()) {
-        /* -------------------------- */
-        /* Cyclic Tridiagonal Section */
-        /* i_r % 2 == 1 */
-        const double h1 = grid.radialSpacing(i_r - 1);
-        const double h2 = grid.radialSpacing(i_r);
-        const double k1 = grid.angularSpacing(i_theta - 1);
-        const double k2 = grid.angularSpacing(i_theta);
-
-        const double coeff1 = 0.5 * (k1 + k2) / h1;
-        const double coeff2 = 0.5 * (k1 + k2) / h2;
-        const double coeff3 = 0.5 * (h1 + h2) / k1;
-        const double coeff4 = 0.5 * (h1 + h2) / k2;
-
-        const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
-        const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
-
-        const int bottom_left  = grid.index(i_r - 1, i_theta_M1);
-        const int left         = grid.index(i_r - 1, i_theta);
-        const int top_left     = grid.index(i_r - 1, i_theta_P1);
-        const int bottom       = grid.index(i_r, i_theta_M1);
-        const int center       = grid.index(i_r, i_theta);
-        const int top          = grid.index(i_r, i_theta_P1);
-        const int bottom_right = grid.index(i_r + 1, i_theta_M1);
-        const int right        = grid.index(i_r + 1, i_theta);
-        const int top_right    = grid.index(i_r + 1, i_theta_P1);
+    KOKKOS_ASSERT(0 < i_r && i_r < grid.numberSmootherCircles());
 
-        if (i_r & 1) {
-            /* i_r % 2 == 1 and i_theta % 2 == 1 */
-            /* | x | o | x | */
-            /* |   |   |   | */
-            /* | o | O | o | */
-            /* |   |   |   | */
-            /* | x | o | x | */
-            /* or */
-            /* i_r % 2 == 1 and i_theta % 2 == 0 */
-            /* | o | o | o | */
-            /* |   |   |   | */
-            /* | x | O | x | */
-            /* |   |   |   | */
-            /* | o | o | o | */
-            result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */
-                                            - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */
-
-                                            - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
-                                            + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
-                                            + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */
-                                            - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */
-                                           );
-        }
-        else {
-            if (i_theta & 1) {
-                /* i_r % 2 == 0 and i_theta % 2 == 1 */
-                /* | o | x | o | */
-                /* |   |   |   | */
-                /* | o | O | o | */
-                /* |   |   |   | */
-                /* | o | x | o | */
-                /* Fill result(i,j) */
-                result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */
-                                                - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */
-                                                - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
-                                                - coeff4 * (att[center] + att[top]) * x[top] /* Top */
-
-                                                - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
-                                                + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
-                                                + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */
-                                                - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */
-                                               );
-            }
-            else {
-                /* i_r % 2 == 0 and i_theta % 2 == 0 */
-                /* | o | o | o | */
-                /* |   |   |   | */
-                /* | o | X | o | */
-                /* |   |   |   | */
-                /* | o | o | o | */
-                result[center] = x[center];
-            }
-        }
-    }
-    /* -------------------- */
-    /* Node on the boundary */
-    /* -------------------- */
-    else if (i_r == 0) {
-        /* ------------------------------------------------ */
-        /* Case 1: Dirichlet boundary on the inner boundary */
-        /* ------------------------------------------------ */
+    if (!(i_r & 1) && !(i_theta & 1)) {
         const int center = grid.index(i_r, i_theta);
-        if (DirBC_Interior) {
-            if (i_theta & 1) {
-                /* i_theta % 2 == 1 */
-                /* || x | o | x | */
-                /* ||   |   |   | */
-                /* || O | o | o | */
-                /* ||   |   |   | */
-                /* || x | o | x | */
-                result[center] = rhs[center];
-            }
-            else {
-                /* i_theta % 2 == 0 */
-                /* || o | o | o | */
-                /* ||   |   |   | */
-                /* || X | o | x | */
-                /* ||   |   |   | */
-                /* || o | o | o | */
-                result[center] = x[center];
-            }
-        }
-        else {
-            /* ------------------------------------------------------------- */
-            /* Case 2: Across origin discretization on the interior boundary */
-            /* ------------------------------------------------------------- */
-            // h1 gets replaced with 2 * R0.
-            // (i_r-1,i_theta) gets replaced with (i_r, i_theta + (grid.ntheta()/2)).
-            // Some more adjustments from the changing the 9-point stencil to the artifical 7-point stencil.
-            const double h1 = 2.0 * grid.radius(0);
-            const double h2 = grid.radialSpacing(i_r);
-            const double k1 = grid.angularSpacing(i_theta - 1);
-            const double k2 = grid.angularSpacing(i_theta);
-
-            const double coeff2 = 0.5 * (k1 + k2) / h2;
-            const double coeff3 = 0.5 * (h1 + h2) / k1;
-            const double coeff4 = 0.5 * (h1 + h2) / k2;
-
-            const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
-            const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
-            //const int i_theta_Across = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2);
-
-            //const int left         = grid.index(i_r, i_theta_Across);
-            const int bottom       = grid.index(i_r, i_theta_M1);
-            const int center       = grid.index(i_r, i_theta);
-            const int top          = grid.index(i_r, i_theta_P1);
-            const int bottom_right = grid.index(i_r + 1, i_theta_M1);
-            const int right        = grid.index(i_r + 1, i_theta);
-            const int top_right    = grid.index(i_r + 1, i_theta_P1);
-
-            if (i_theta & 1) {
-                /* i_theta % 2 == 1 */
-                /* -| x | o | x | */
-                /* -|   |   |   | */
-                /* -| O | o | o | */
-                /* -|   |   |   | */
-                /* -| x | o | x | */
-                result[center] =
-                    rhs[center] -
-                    (-coeff2 * (arr[center] + arr[right]) * x[right] /* Right */
-                     - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
-                     - coeff4 * (att[center] + att[top]) * x[top] /* Top */
-
-                     /* - 0.25 * (art[left] + art[bottom]) * x[bottom_left] // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
-                     + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
-
-                     /* + 0.25 * (art[left] + art[top]) * x[top_left] // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
-                     - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */
-                    );
-            }
-            else {
-                /* i_theta % 2 == 0 */
-                /* -| o | o | o | */
-                /* -|   |   |   | */
-                /* -| X | o | x | */
-                /* -|   |   |   | */
-                /* -| o | o | o | */
-                result[center] = x[center];
-            }
-        }
+        result[center]   = x[center];
+        return;
+    }
+
+    const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
+    const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
+
+    const int bottom_left  = grid.index(i_r - 1, i_theta_M1);
+    const int left         = grid.index(i_r - 1, i_theta);
+    const int top_left     = grid.index(i_r - 1, i_theta_P1);
+    const int bottom       = grid.index(i_r, i_theta_M1);
+    const int center       = grid.index(i_r, i_theta);
+    const int top          = grid.index(i_r, i_theta_P1);
+    const int bottom_right = grid.index(i_r + 1, i_theta_M1);
+    const int right        = grid.index(i_r + 1, i_theta);
+    const int top_right    = grid.index(i_r + 1, i_theta_P1);
+
+    const double h1 = grid.radialSpacing(i_r - 1);
+    const double h2 = grid.radialSpacing(i_r);
+    const double k1 = grid.angularSpacing(i_theta_M1);
+    const double k2 = grid.angularSpacing(i_theta);
+
+    const double coeff1 = 0.5 * (k1 + k2) / h1;
+    const double coeff2 = 0.5 * (k1 + k2) / h2;
+    const double coeff3 = 0.5 * (h1 + h2) / k1;
+    const double coeff4 = 0.5 * (h1 + h2) / k2;
+
+    result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */
+                                    - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */
+
+                                    - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
+                                    + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
+                                    + (art[left] + art[top]) * x[top_left] /* Top Left */
+                                    - (art[right] + art[top]) * x[top_right] /* Top Right */
+                                   );
+
+    if (!(i_r & 1) && (i_theta & 1)) {
+        /* | o | x | o | */
+        /* |   |   |   | */
+        /* | o | O | o | */
+        /* |   |   |   | */
+        /* | o | x | o | */
+        result[center] -= (-coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
+                           - coeff4 * (att[center] + att[top]) * x[top] /* Top */
+        );
     }
 }
 
-static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(const int i_r, const int i_theta, const PolarGrid& grid,
-                                                               const bool DirBC_Interior, ConstVector<double>& x,
-                                                               ConstVector<double>& rhs, Vector<double>& result,
-                                                               ConstVector<double>& arr, ConstVector<double>& att,
-                                                               ConstVector<double>& art)
+static KOKKOS_INLINE_FUNCTION void
+nodeApplyAscOrthoCircleTakeBoundary(const int i_r, const int i_theta, const PolarGrid& grid, const bool DirBC_Interior,
+                                    ConstVector<double>& x, ConstVector<double>& rhs, Vector<double>& result,
+                                    ConstVector<double>& arr, ConstVector<double>& att, ConstVector<double>& art)
 {
-    assert(i_r >= grid.numberSmootherCircles() - 1 && i_r < grid.nr());
+    KOKKOS_ASSERT(i_r == 0);
 
-    /* -------------------- */
-    /* Node in the interior */
-    /* -------------------- */
-    if (i_r > grid.numberSmootherCircles() && i_r < grid.nr() - 2) {
-        const double h1 = grid.radialSpacing(i_r - 1);
+    if (!(i_theta & 1)) {
+        const int center = grid.index(i_r, i_theta);
+        result[center]   = x[center];
+        return;
+    }
+
+    if (DirBC_Interior) {
+        const int center = grid.index(i_r, i_theta);
+        result[center]   = rhs[center];
+    }
+    else {
+        const double h1 = 2.0 * grid.radius(0);
         const double h2 = grid.radialSpacing(i_r);
         const double k1 = grid.angularSpacing(i_theta - 1);
         const double k2 = grid.angularSpacing(i_theta);
 
-        const double coeff1 = 0.5 * (k1 + k2) / h1;
         const double coeff2 = 0.5 * (k1 + k2) / h2;
         const double coeff3 = 0.5 * (h1 + h2) / k1;
         const double coeff4 = 0.5 * (h1 + h2) / k2;
@@ -209,9 +90,6 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(const int i_r, co
         const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
         const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
 
-        const int bottom_left  = grid.index(i_r - 1, i_theta_M1);
-        const int left         = grid.index(i_r - 1, i_theta);
-        const int top_left     = grid.index(i_r - 1, i_theta_P1);
         const int bottom       = grid.index(i_r, i_theta_M1);
         const int center       = grid.index(i_r, i_theta);
         const int top          = grid.index(i_r, i_theta_P1);
@@ -219,146 +97,139 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(const int i_r, co
         const int right        = grid.index(i_r + 1, i_theta);
         const int top_right    = grid.index(i_r + 1, i_theta_P1);
 
-        if (i_theta & 1) {
-            /* i_theta % 2 == 1 and i_r % 2 == 1 */
-            /* ---------- */
-            /* x   o   x  */
-            /* ---------- */
-            /* o   O   o  */
-            /* ---------- */
-            /* x   o   x  */
-            /* ---------- */
-            /* or */
-            /* i_theta % 2 == 1 and i_r % 2 == 0 */
-            /* ---------- */
-            /* o   x   o  */
-            /* ---------- */
-            /* o   O   o  */
-            /* ---------- */
-            /* o   x   o  */
-            /* ---------- */
-            result[center] = rhs[center] - (-coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
-                                            - coeff4 * (att[center] + att[top]) * x[top] /* Top */
+        result[center] = rhs[center] - (-coeff2 * (arr[center] + arr[right]) * x[right] /* Right */
+                                        - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
+                                        - coeff4 * (att[center] + att[top]) * x[top] /* Top */
 
-                                            - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
-                                            + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
-                                            + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */
-                                            - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */
-                                           );
-        }
-        else {
-            if (i_r & 1) {
-                /* i_theta % 2 == 0 and i_r % 2 == 1 */
-                /* ---------- */
-                /* o   o   o  */
-                /* ---------- */
-                /* x   O   x  */
-                /* ---------- */
-                /* o   o   o  */
-                /* ---------- */
-                result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */
-                                                - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */
-                                                - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
-                                                - coeff4 * (att[center] + att[top]) * x[top] /* Top */
-
-                                                - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
-                                                + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
-                                                + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */
-                                                - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */
-                                               );
-            }
-            else {
-                /* i_theta % 2 == 0 and i_r % 2 == 0 */
-                /* ---------- */
-                /* o   o   o  */
-                /* ---------- */
-                /* o   X   o  */
-                /* ---------- */
-                /* o   o   o  */
-                /* ---------- */
-                result[center] = x[center];
-            }
-        }
+                                        + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
+                                        - (art[right] + art[top]) * x[top_right] /* Top Right */
+                                       );
     }
-    else if (i_r == grid.numberSmootherCircles()) {
-        const double h1 = grid.radialSpacing(i_r - 1);
-        const double h2 = grid.radialSpacing(i_r);
-        const double k1 = grid.angularSpacing(i_theta - 1);
-        const double k2 = grid.angularSpacing(i_theta);
+}
 
-        const double coeff1 = 0.5 * (k1 + k2) / h1;
-        const double coeff2 = 0.5 * (k1 + k2) / h2;
-        const double coeff3 = 0.5 * (h1 + h2) / k1;
-        const double coeff4 = 0.5 * (h1 + h2) / k2;
+static KOKKOS_INLINE_FUNCTION void
+nodeApplyAscOrthoRadialTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, ConstVector<double>& x,
+                                    ConstVector<double>& rhs, Vector<double>& result, ConstVector<double>& arr,
+                                    ConstVector<double>& att, ConstVector<double>& art)
+{
+    KOKKOS_ASSERT(grid.numberSmootherCircles() < i_r && i_r < grid.nr() - 2);
 
-        const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
-        const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
+    if (!(i_r & 1) && !(i_theta & 1)) {
+        const int center = grid.index(i_r, i_theta);
+        result[center]   = x[center];
+        return;
+    }
 
-        const int bottom_left  = grid.index(i_r - 1, i_theta_M1);
-        const int left         = grid.index(i_r - 1, i_theta);
-        const int top_left     = grid.index(i_r - 1, i_theta_P1);
-        const int bottom       = grid.index(i_r, i_theta_M1);
-        const int center       = grid.index(i_r, i_theta);
-        const int top          = grid.index(i_r, i_theta_P1);
-        const int bottom_right = grid.index(i_r + 1, i_theta_M1);
-        const int right        = grid.index(i_r + 1, i_theta);
-        const int top_right    = grid.index(i_r + 1, i_theta_P1);
+    const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
+    const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
+
+    const int bottom_left  = grid.index(i_r - 1, i_theta_M1);
+    const int bottom       = grid.index(i_r, i_theta_M1);
+    const int bottom_right = grid.index(i_r + 1, i_theta_M1);
+    const int left         = grid.index(i_r - 1, i_theta);
+    const int center       = grid.index(i_r, i_theta);
+    const int right        = grid.index(i_r + 1, i_theta);
+    const int top_left     = grid.index(i_r - 1, i_theta_P1);
+    const int top          = grid.index(i_r, i_theta_P1);
+    const int top_right    = grid.index(i_r + 1, i_theta_P1);
+
+    const double h1 = grid.radialSpacing(i_r - 1);
+    const double h2 = grid.radialSpacing(i_r);
+    const double k1 = grid.angularSpacing(i_theta - 1);
+    const double k2 = grid.angularSpacing(i_theta);
+
+    const double coeff1 = 0.5 * (k1 + k2) / h1;
+    const double coeff2 = 0.5 * (k1 + k2) / h2;
+    const double coeff3 = 0.5 * (h1 + h2) / k1;
+    const double coeff4 = 0.5 * (h1 + h2) / k2;
+
+    result[center] = rhs[center] - (-coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
+                                    - coeff4 * (att[center] + att[top]) * x[top] /* Top */
+
+                                    - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
+                                    + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
+                                    + (art[left] + art[top]) * x[top_left] /* Top Left */
+                                    - (art[right] + art[top]) * x[top_right] /* Top Right */
+                                   );
+
+    if (!(i_theta & 1) && (i_r & 1)) {
+        /* ---------- */
+        /* o   o   o  */
+        /* ---------- */
+        /* x   O   x  */
+        /* ---------- */
+        /* o   o   o  */
+        /* ---------- */
+        result[center] -= (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */
+                           - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */
+        );
+    }
+}
 
-        if (i_theta & 1) {
-            /* i_theta % 2 == 1 and i_r % 2 == 1 */
-            /* | x | o | x || o   x   o   x  */
-            /* |   |   |   || -------------- */
-            /* | o | o | o || O   o   o   o  */
-            /* |   |   |   || -------------- */
-            /* | x | o | x || o   x   o   x  */
-            /* or */
-            /* i_theta % 2 == 1 and i_r % 2 == 0 */
-            /* | o | x | o || x   o   x   o  */
-            /* |   |   |   || -------------- */
-            /* | o | o | o || O   o   o   o  */
-            /* |   |   |   || -------------- */
-            /* | o | x | o || x   o   x   o  */
-            result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */
-                                            - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
-                                            - coeff4 * (att[center] + att[top]) * x[top] /* Top */
+static KOKKOS_INLINE_FUNCTION void
+nodeApplyAscOrthoRadialTakeInnerBoundary(const int i_r, const int i_theta, const PolarGrid& grid,
+                                         ConstVector<double>& x, ConstVector<double>& rhs, Vector<double>& result,
+                                         ConstVector<double>& arr, ConstVector<double>& att, ConstVector<double>& art)
+{
+    KOKKOS_ASSERT(i_r == grid.numberSmootherCircles());
 
-                                            - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
-                                            + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
-                                            + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */
-                                            - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */
-                                           );
-        }
-        else {
-            if (i_r & 1) {
-                /* i_theta % 2 == 0 and i_r % 2 == 1 */
-                /* | o | o | o || o   o   o   o  */
-                /* |   |   |   || -------------- */
-                /* | x | o | x || O   x   o   x  */
-                /* |   |   |   || -------------- */
-                /* | o | o | o || o   o   o   o  */
-                result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */
-                                                - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */
-                                                - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
-                                                - coeff4 * (att[center] + att[top]) * x[top] /* Top */
-
-                                                - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
-                                                + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
-                                                + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */
-                                                - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */
-                                               );
-            }
-            else {
-                /* i_theta % 2 == 0 and i_r % 2 == 0 */
-                /* | o | o | o || o   o   o   o  */
-                /* |   |   |   || -------------- */
-                /* | o | x | o || X   o   x   o  */
-                /* |   |   |   || -------------- */
-                /* | o | o | o || o   o   o   o  */
-                result[center] = x[center];
-            }
-        }
+    if (!(i_r & 1) && !(i_theta & 1)) {
+        const int center = grid.index(i_r, i_theta);
+        result[center]   = x[center];
+        return;
     }
-    else if (i_r == grid.nr() - 2) {
+
+    const double h1 = grid.radialSpacing(i_r - 1);
+    const double h2 = grid.radialSpacing(i_r);
+    const double k1 = grid.angularSpacing(i_theta - 1);
+    const double k2 = grid.angularSpacing(i_theta);
+
+    const double coeff1 = 0.5 * (k1 + k2) / h1;
+    const double coeff2 = 0.5 * (k1 + k2) / h2;
+    const double coeff3 = 0.5 * (h1 + h2) / k1;
+    const double coeff4 = 0.5 * (h1 + h2) / k2;
+
+    const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
+    const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
+
+    const int bottom_left  = grid.index(i_r - 1, i_theta_M1);
+    const int left         = grid.index(i_r - 1, i_theta);
+    const int top_left     = grid.index(i_r - 1, i_theta_P1);
+    const int bottom       = grid.index(i_r, i_theta_M1);
+    const int center       = grid.index(i_r, i_theta);
+    const int top          = grid.index(i_r, i_theta_P1);
+    const int bottom_right = grid.index(i_r + 1, i_theta_M1);
+    const int right        = grid.index(i_r + 1, i_theta);
+    const int top_right    = grid.index(i_r + 1, i_theta_P1);
+
+    result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */
+                                    - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
+                                    - coeff4 * (att[center] + att[top]) * x[top] /* Top */
+
+                                    - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
+                                    + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
+                                    + (art[left] + art[top]) * x[top_left] /* Top Left */
+                                    - (art[right] + art[top]) * x[top_right] /* Top Right */
+                                   );
+
+    if (!(i_theta & 1) && (i_r & 1)) {
+        /* | o | o | o || o   o   o   o  */
+        /* |   |   |   || -------------- */
+        /* | x | o | x || O   x   o   x  */
+        /* |   |   |   || -------------- */
+        /* | o | o | o || o   o   o   o  */
+        result[center] -= (-coeff2 * (arr[center] + arr[right]) * x[right] /* Right */);
+    }
+}
+
+static KOKKOS_INLINE_FUNCTION void
+nodeApplyAscOrthoRadialTakeOuterBoundary(const int i_r, const int i_theta, const PolarGrid& grid,
+                                         ConstVector<double>& x, ConstVector<double>& rhs, Vector<double>& result,
+                                         ConstVector<double>& arr, ConstVector<double>& att, ConstVector<double>& art)
+{
+    KOKKOS_ASSERT(i_r == grid.nr() - 2 || i_r == grid.nr() - 1);
+
+    if (i_r == grid.nr() - 2) {
         KOKKOS_ASSERT(i_r & 1);
 
         const double h1 = grid.radialSpacing(i_r - 1);
@@ -385,7 +256,6 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(const int i_r, co
         const int top_right    = grid.index(i_r + 1, i_theta_P1);
 
         if (i_theta & 1) {
-            /* i_theta % 2 == 1 */
             /* ---------------|| */
             /* o   x   o   x  || */
             /* ---------------|| */
@@ -402,10 +272,10 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(const int i_r, co
                                - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
                                - coeff4 * (att[center] + att[top]) * x[top] /* Top */
 
-                               - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
-                               + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
-                               + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */
-                               - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */
+                               - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
+                               + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
+                               + (art[left] + art[top]) * x[top_left] /* Top Left */
+                               - (art[right] + art[top]) * x[top_right] /* Top Right */
                               );
         }
         else {
@@ -421,10 +291,10 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(const int i_r, co
                                             - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
                                             - coeff4 * (att[center] + att[top]) * x[top] /* Top */
 
-                                            - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
-                                            + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
-                                            + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */
-                                            - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */
+                                            - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
+                                            + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
+                                            + (art[left] + art[top]) * x[top_left] /* Top Left */
+                                            - (art[right] + art[top]) * x[top_right] /* Top Right */
                                            );
         }
     }
@@ -434,7 +304,6 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(const int i_r, co
         const int center = grid.index(i_r, i_theta);
 
         if (i_theta & 1) {
-            /* i_theta % 2 == 1 */
             /* -----------|| */
             /* x   o   x  || */
             /* -----------|| */
@@ -464,7 +333,8 @@ void ExtrapolatedSmootherTake<LevelCacheType>::applyAscOrthoBlackCircleSection(C
                                                                                ConstVector<double> rhs,
                                                                                Vector<double> temp)
 {
-    using extrapolated_smoother_take::nodeApplyAscOrthoCircleTake;
+    using extrapolated_smoother_take::nodeApplyAscOrthoCircleTakeBoundary;
+    using extrapolated_smoother_take::nodeApplyAscOrthoCircleTakeInterior;
 
     const PolarGrid& grid             = ExtrapolatedSmoother<LevelCacheType>::grid_;
     const LevelCacheType& level_cache = ExtrapolatedSmoother<LevelCacheType>::level_cache_;
@@ -480,17 +350,44 @@ void ExtrapolatedSmootherTake<LevelCacheType>::applyAscOrthoBlackCircleSection(C
     const int start_black_circles = (grid.numberSmootherCircles() % 2 == 0) ? 1 : 0;
     const int num_black_circles   = (grid.numberSmootherCircles() - start_black_circles + 1) / 2;
 
-    Kokkos::parallel_for(
-        "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Circular)",
-        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
-            {0, 0}, // Starting point of the index space
-            {num_black_circles, grid.ntheta()} // Ending point of the index space
-            ),
-        // Kokkos lambda function to execute for each point in the index space
-        KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
-            int i_r = start_black_circles + circle_task * 2;
-            nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art);
-        });
+    if (start_black_circles == 0) {
+        Kokkos::parallel_for(
+            "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Circular)",
+            Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+                {0, 0}, // Starting point of the index space
+                {1, grid.ntheta()} // Ending point of the index space
+                ),
+            // Kokkos lambda function to execute for each point in the index space
+            KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
+                int i_r = start_black_circles + circle_task * 2;
+                nodeApplyAscOrthoCircleTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art);
+            });
+
+        Kokkos::parallel_for(
+            "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Circular)",
+            Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+                {1, 0}, // Starting point of the index space
+                {num_black_circles, grid.ntheta()} // Ending point of the index space
+                ),
+            // Kokkos lambda function to execute for each point in the index space
+            KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
+                int i_r = start_black_circles + circle_task * 2;
+                nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+            });
+    }
+    else {
+        Kokkos::parallel_for(
+            "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Circular)",
+            Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+                {0, 0}, // Starting point of the index space
+                {num_black_circles, grid.ntheta()} // Ending point of the index space
+                ),
+            // Kokkos lambda function to execute for each point in the index space
+            KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
+                int i_r = start_black_circles + circle_task * 2;
+                nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+            });
+    }
 
     Kokkos::fence();
 }
@@ -500,7 +397,8 @@ void ExtrapolatedSmootherTake<LevelCacheType>::applyAscOrthoWhiteCircleSection(C
                                                                                ConstVector<double> rhs,
                                                                                Vector<double> temp)
 {
-    using extrapolated_smoother_take::nodeApplyAscOrthoCircleTake;
+    using extrapolated_smoother_take::nodeApplyAscOrthoCircleTakeBoundary;
+    using extrapolated_smoother_take::nodeApplyAscOrthoCircleTakeInterior;
 
     const PolarGrid& grid             = ExtrapolatedSmoother<LevelCacheType>::grid_;
     const LevelCacheType& level_cache = ExtrapolatedSmoother<LevelCacheType>::level_cache_;
@@ -516,17 +414,44 @@ void ExtrapolatedSmootherTake<LevelCacheType>::applyAscOrthoWhiteCircleSection(C
     const int start_white_circles = (grid.numberSmootherCircles() % 2 == 0) ? 0 : 1;
     const int num_white_circles   = (grid.numberSmootherCircles() - start_white_circles + 1) / 2;
 
-    Kokkos::parallel_for(
-        "ExtrapolatedSmootherTake: ApplyAscOrtho (White Circular)",
-        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
-            {0, 0}, // Starting point of the index space
-            {num_white_circles, grid.ntheta()} // Ending point of the index space
-            ),
-        // Kokkos lambda function to execute for each point in the index space
-        KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
-            const int i_r = start_white_circles + circle_task * 2;
-            nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art);
-        });
+    if (start_white_circles == 0) {
+        Kokkos::parallel_for(
+            "ExtrapolatedSmootherTake: ApplyAscOrtho (White Circular)",
+            Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+                {0, 0}, // Starting point of the index space
+                {1, grid.ntheta()} // Ending point of the index space
+                ),
+            // Kokkos lambda function to execute for each point in the index space
+            KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
+                const int i_r = start_white_circles + circle_task * 2;
+                nodeApplyAscOrthoCircleTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art);
+            });
+
+        Kokkos::parallel_for(
+            "ExtrapolatedSmootherTake: ApplyAscOrtho (White Circular)",
+            Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+                {1, 0}, // Starting point of the index space
+                {num_white_circles, grid.ntheta()} // Ending point of the index space
+                ),
+            // Kokkos lambda function to execute for each point in the index space
+            KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
+                const int i_r = start_white_circles + circle_task * 2;
+                nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+            });
+    }
+    else {
+        Kokkos::parallel_for(
+            "ExtrapolatedSmootherTake: ApplyAscOrtho (White Circular)",
+            Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+                {0, 0}, // Starting point of the index space
+                {num_white_circles, grid.ntheta()} // Ending point of the index space
+                ),
+            // Kokkos lambda function to execute for each point in the index space
+            KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
+                const int i_r = start_white_circles + circle_task * 2;
+                nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+            });
+    }
 
     Kokkos::fence();
 }
@@ -536,11 +461,12 @@ void ExtrapolatedSmootherTake<LevelCacheType>::applyAscOrthoBlackRadialSection(C
                                                                                ConstVector<double> rhs,
                                                                                Vector<double> temp)
 {
-    using extrapolated_smoother_take::nodeApplyAscOrthoRadialTake;
+    using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeInnerBoundary;
+    using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeInterior;
+    using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeOuterBoundary;
 
     const PolarGrid& grid             = ExtrapolatedSmoother<LevelCacheType>::grid_;
     const LevelCacheType& level_cache = ExtrapolatedSmoother<LevelCacheType>::level_cache_;
-    const bool DirBC_Interior         = ExtrapolatedSmoother<LevelCacheType>::DirBC_Interior_;
 
     assert(level_cache.cacheDomainGeometry());
 
@@ -556,12 +482,36 @@ void ExtrapolatedSmootherTake<LevelCacheType>::applyAscOrthoBlackRadialSection(C
         "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Radial)",
         Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
             {0, grid.numberSmootherCircles()}, // Starting point of the index space
+            {num_black_radial_lines, grid.numberSmootherCircles() + 1} // Ending point of the index space
+            ),
+        // Kokkos lambda function to execute for each point in the index space
+        KOKKOS_LAMBDA(const int radial_task, const int i_r) {
+            const int i_theta = start_black_radials + radial_task * 2;
+            nodeApplyAscOrthoRadialTakeInnerBoundary(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+        });
+
+    Kokkos::parallel_for(
+        "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Radial)",
+        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+            {0, grid.numberSmootherCircles() + 1}, // Starting point of the index space
+            {num_black_radial_lines, grid.nr() - 2} // Ending point of the index space
+            ),
+        // Kokkos lambda function to execute for each point in the index space
+        KOKKOS_LAMBDA(const int radial_task, const int i_r) {
+            const int i_theta = start_black_radials + radial_task * 2;
+            nodeApplyAscOrthoRadialTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+        });
+
+    Kokkos::parallel_for(
+        "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Radial)",
+        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+            {0, grid.nr() - 2}, // Starting point of the index space
             {num_black_radial_lines, grid.nr()} // Ending point of the index space
             ),
         // Kokkos lambda function to execute for each point in the index space
         KOKKOS_LAMBDA(const int radial_task, const int i_r) {
             const int i_theta = start_black_radials + radial_task * 2;
-            nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art);
+            nodeApplyAscOrthoRadialTakeOuterBoundary(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
         });
 
     Kokkos::fence();
@@ -572,11 +522,12 @@ void ExtrapolatedSmootherTake<LevelCacheType>::applyAscOrthoWhiteRadialSection(C
                                                                                ConstVector<double> rhs,
                                                                                Vector<double> temp)
 {
-    using extrapolated_smoother_take::nodeApplyAscOrthoRadialTake;
+    using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeInnerBoundary;
+    using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeInterior;
+    using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeOuterBoundary;
 
     const PolarGrid& grid             = ExtrapolatedSmoother<LevelCacheType>::grid_;
     const LevelCacheType& level_cache = ExtrapolatedSmoother<LevelCacheType>::level_cache_;
-    const bool DirBC_Interior         = ExtrapolatedSmoother<LevelCacheType>::DirBC_Interior_;
 
     assert(level_cache.cacheDomainGeometry());
 
@@ -592,12 +543,36 @@ void ExtrapolatedSmootherTake<LevelCacheType>::applyAscOrthoWhiteRadialSection(C
         "ExtrapolatedSmootherTake: ApplyAscOrtho (White Radial)",
         Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
             {0, grid.numberSmootherCircles()}, // Starting point of the index space
+            {num_white_radial_lines, grid.numberSmootherCircles() + 1} // Ending point of the index space
+            ),
+        // Kokkos lambda function to execute for each point in the index space
+        KOKKOS_LAMBDA(const int radial_task, const int i_r) {
+            const int i_theta = start_white_radials + radial_task * 2;
+            nodeApplyAscOrthoRadialTakeInnerBoundary(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+        });
+
+    Kokkos::parallel_for(
+        "ExtrapolatedSmootherTake: ApplyAscOrtho (White Radial)",
+        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+            {0, grid.numberSmootherCircles() + 1}, // Starting point of the index space
+            {num_white_radial_lines, grid.nr() - 2} // Ending point of the index space
+            ),
+        // Kokkos lambda function to execute for each point in the index space
+        KOKKOS_LAMBDA(const int radial_task, const int i_r) {
+            const int i_theta = start_white_radials + radial_task * 2;
+            nodeApplyAscOrthoRadialTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+        });
+
+    Kokkos::parallel_for(
+        "ExtrapolatedSmootherTake: ApplyAscOrtho (White Radial)",
+        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+            {0, grid.nr() - 2}, // Starting point of the index space
             {num_white_radial_lines, grid.nr()} // Ending point of the index space
             ),
         // Kokkos lambda function to execute for each point in the index space
         KOKKOS_LAMBDA(const int radial_task, const int i_r) {
             const int i_theta = start_white_radials + radial_task * 2;
-            nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art);
+            nodeApplyAscOrthoRadialTakeOuterBoundary(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
         });
 
     Kokkos::fence();
diff --git a/include/Residual/ResidualGive/applyAGive.inl b/include/Residual/ResidualGive/applyAGive.inl
index 7370f37e..d016a5e0 100644
--- a/include/Residual/ResidualGive/applyAGive.inl
+++ b/include/Residual/ResidualGive/applyAGive.inl
@@ -53,26 +53,26 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_give(int i_r, int i_theta, const
         if (i_r > 1 || (i_r == 1 && !DirBC_Interior)) { /* Don't give to the inner dirichlet boundary! */
             result[left] += (-coeff1 * arr * x[center] /* Right */
                              + coeff1 * arr * x[left] /* Center: (Right) */
-                             - 0.25 * art * x[top] /* Top Right */
-                             + 0.25 * art * x[bottom]); /* Bottom Right */
+                             - art * x[top] /* Top Right */
+                             + art * x[bottom]); /* Bottom Right */
         }
         /* Fill result(i+1,j) */
         if (i_r < grid.nr() - 2) { /* Don't give to the outer dirichlet boundary! */
             result[right] += (-coeff2 * arr * x[center] /* Left */
                               + coeff2 * arr * x[right] /* Center: (Left) */
-                              + 0.25 * art * x[top] /* Top Left */
-                              - 0.25 * art * x[bottom]); /* Bottom Left */
+                              + art * x[top] /* Top Left */
+                              - art * x[bottom]); /* Bottom Left */
         }
         /* Fill result(i,j-1) */
         result[bottom] += (-coeff3 * att * x[center] /* Top */
                            + coeff3 * att * x[bottom] /* Center: (Top) */
-                           - 0.25 * art * x[right] /* Top Right */
-                           + 0.25 * art * x[left]); /* Top Left */
+                           - art * x[right] /* Top Right */
+                           + art * x[left]); /* Top Left */
         /* Fill result(i,j+1) */
         result[top] += (-coeff4 * att * x[center] /* Bottom */
                         + coeff4 * att * x[top] /* Center: (Bottom) */
-                        + 0.25 * art * x[right] /* Bottom Right */
-                        - 0.25 * art * x[left]); /* Bottom Left */
+                        + art * x[right] /* Bottom Right */
+                        - art * x[left]); /* Bottom Left */
     }
     /* -------------------------- */
     /* Node on the inner boundary */
@@ -102,8 +102,8 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_give(int i_r, int i_theta, const
             /* Fill result(i+1,j) */
             result[right] += (-coeff2 * arr * x[center] /* Left */
                               + coeff2 * arr * x[right] /* Center: (Left) */
-                              + 0.25 * art * x[top] /* Top Left */
-                              - 0.25 * art * x[bottom]); /* Bottom Left */
+                              + art * x[top] /* Top Left */
+                              - art * x[bottom]); /* Bottom Left */
         }
         else {
             /* ------------------------------------------------------------- */
@@ -144,23 +144,23 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_give(int i_r, int i_theta, const
             /* From view the view of the across origin node, the directions are roatated by 180 degrees in the stencil! */
             result[left] += (-coeff1 * arr * x[center] /* Right -> Left */
                              + coeff1 * arr * x[left]); /* Center: (Right) -> Center: (Left)*/
-            /* + 0.25 * art * x[top]; // Top Right -> Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
-            /* - 0.25 * art * x[bottom]; // Bottom Right -> Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
+            /* + art * x[top]; // Top Right -> Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
+            /* - art * x[bottom]; // Bottom Right -> Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
             /* Fill result(i+1,j) */
             result[right] += (-coeff2 * arr * x[center] /* Left */
                               + coeff2 * arr * x[right] /* Center: (Left) */
-                              + 0.25 * art * x[top] /* Top Left */
-                              - 0.25 * art * x[bottom]); /* Bottom Left */
+                              + art * x[top] /* Top Left */
+                              - art * x[bottom]); /* Bottom Left */
             /* Fill result(i,j-1) */
             result[bottom] += (-coeff3 * att * x[center] /* Top */
                                + coeff3 * att * x[bottom] /* Center: (Top) */
-                               - 0.25 * art * x[right]); /* Top Right */
-            /* + 0.25 * art * x[left]; // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
+                               - art * x[right]); /* Top Right */
+            /* + art * x[left]; // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
             /* Fill result(i,j+1) */
             result[top] += (-coeff4 * att * x[center] /* Bottom */
                             + coeff4 * att * x[top] /* Center: (Bottom) */
-                            + 0.25 * art * x[right]); /* Bottom Right */
-            /* - 0.25 * art * x[left]; // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
+                            + art * x[right]); /* Bottom Right */
+            /* - art * x[left]; // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
         }
     }
     /* ----------------------------- */
@@ -187,8 +187,8 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_give(int i_r, int i_theta, const
         /* Fill result(i-1,j) */
         result[left] += (-coeff1 * arr * x[center] /* Right */
                          + coeff1 * arr * x[left] /* Center: (Right) */
-                         - 0.25 * art * x[top] /* Top Right */
-                         + 0.25 * art * x[bottom]); /* Bottom Right */
+                         - art * x[top] /* Top Right */
+                         + art * x[bottom]); /* Bottom Right */
     }
 }
 
diff --git a/include/Residual/ResidualGive/residualGive.h b/include/Residual/ResidualGive/residualGive.h
index cb51ad8f..595c3b9b 100644
--- a/include/Residual/ResidualGive/residualGive.h
+++ b/include/Residual/ResidualGive/residualGive.h
@@ -10,7 +10,6 @@ class ResidualGive : public Residual<LevelCacheType>
 {
 public:
     explicit ResidualGive(const PolarGrid& grid, const LevelCacheType& level_cache, const bool DirBC_Interior);
-    ~ResidualGive() override = default;
 
     void applySystemOperator(Vector<double> result, ConstVector<double> x) const final;
     void computeResidual(Vector<double> result, ConstVector<double> rhs, ConstVector<double> x) const final;
diff --git a/include/Residual/ResidualTake/applyATake.inl b/include/Residual/ResidualTake/applyATake.inl
index 9d39f520..6eb47a29 100644
--- a/include/Residual/ResidualTake/applyATake.inl
+++ b/include/Residual/ResidualTake/applyATake.inl
@@ -3,23 +3,29 @@
 namespace residual_take
 {
 
-static KOKKOS_INLINE_FUNCTION void node_apply_a_take(const int i_r, const int i_theta, const PolarGrid& grid,
-                                                     bool DirBC_Interior, Vector<double>& result,
-                                                     ConstVector<double>& x, ConstVector<double>& arr,
-                                                     ConstVector<double>& att, ConstVector<double>& art,
-                                                     ConstVector<double>& detDF, ConstVector<double>& coeff_beta)
+static KOKKOS_INLINE_FUNCTION void
+applySystemMatrixTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, Vector<double>& result,
+                              ConstVector<double>& x, ConstVector<double>& arr, ConstVector<double>& att,
+                              ConstVector<double>& art, ConstVector<double>& detDF, ConstVector<double>& coeff_beta)
 {
-    const int center = grid.index(i_r, i_theta);
+    KOKKOS_ASSERT(0 < i_r && i_r < grid.nr() - 1);
 
-    if ((i_r == 0 && DirBC_Interior) || (i_r == grid.nr() - 1)) {
-        result[center] = x[center];
-        return;
-    }
+    const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
+    const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
 
-    // Across origin: h1 gets replaced with 2 * R0.
-    const double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1);
+    const int bottom_left  = grid.index(i_r - 1, i_theta_M1);
+    const int bottom       = grid.index(i_r, i_theta_M1);
+    const int bottom_right = grid.index(i_r + 1, i_theta_M1);
+    const int left         = grid.index(i_r - 1, i_theta);
+    const int center       = grid.index(i_r, i_theta);
+    const int right        = grid.index(i_r + 1, i_theta);
+    const int top_left     = grid.index(i_r - 1, i_theta_P1);
+    const int top          = grid.index(i_r, i_theta_P1);
+    const int top_right    = grid.index(i_r + 1, i_theta_P1);
+
+    const double h1 = grid.radialSpacing(i_r - 1);
     const double h2 = grid.radialSpacing(i_r);
-    const double k1 = grid.angularSpacing(i_theta - 1);
+    const double k1 = grid.angularSpacing(i_theta_M1);
     const double k2 = grid.angularSpacing(i_theta);
 
     const double coeff1 = 0.5 * (k1 + k2) / h1;
@@ -28,38 +34,69 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_take(const int i_r, const int i_
     const double coeff4 = 0.5 * (h1 + h2) / k2;
     const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2);
 
-    const int i_theta_M1     = grid.wrapThetaIndex(i_theta - 1);
-    const int i_theta_P1     = grid.wrapThetaIndex(i_theta + 1);
-    const int i_theta_Across = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2);
-
-    const int left   = (i_r == 0) ? grid.index(i_r, i_theta_Across) : grid.index(i_r - 1, i_theta);
-    const int right  = grid.index(i_r + 1, i_theta);
-    const int bottom = grid.index(i_r, i_theta_M1);
-    const int top    = grid.index(i_r, i_theta_P1);
+    result[center] = (+coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) * x[center] /* beta_{i,j} */
 
-    double value = 0.0;
+                      + coeff1 * (arr[center] + arr[left]) * (x[center] - x[left]) /* Center: (Left) - Left */
+                      + coeff2 * (arr[center] + arr[right]) * (x[center] - x[right]) /* Center: (Right) - Right */
+                      + coeff3 * (att[center] + att[bottom]) * (x[center] - x[bottom]) /* Center: (Bottom) - Bottom */
+                      + coeff4 * (att[center] + att[top]) * (x[center] - x[top]) /* Center: (Top) - Top */
 
-    value += coeff5 * coeff_beta[center] * std::fabs(detDF[center]) * x[center]; /* beta_{i,j} */
+                      - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
+                      + (art[left] + art[top]) * x[top_left] /* Top Left */
+                      + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
+                      - (art[right] + art[top]) * x[top_right] /* Top Right */
+    );
+}
 
-    value += coeff1 * (arr[center] + arr[left]) * (x[center] - x[left]); /* Center: (Left) - Left */
-    value += coeff2 * (arr[center] + arr[right]) * (x[center] - x[right]); /* Center: (Right) - Right */
-    value += coeff3 * (att[center] + att[bottom]) * (x[center] - x[bottom]); /* Center: (Bottom) - Bottom */
-    value += coeff4 * (att[center] + att[top]) * (x[center] - x[top]); /* Center: (Top) - Top */
+static KOKKOS_INLINE_FUNCTION void applySystemMatrixTakeBoundary(const int i_r, const int i_theta,
+                                                                 const PolarGrid& grid, bool DirBC_Interior,
+                                                                 Vector<double>& result, ConstVector<double>& x,
+                                                                 ConstVector<double>& arr, ConstVector<double>& att,
+                                                                 ConstVector<double>& art, ConstVector<double>& detDF,
+                                                                 ConstVector<double>& coeff_beta)
+{
+    KOKKOS_ASSERT(i_r == 0 || i_r == grid.nr() - 1);
 
-    // Across origin: Reduce 9-point stencil to the artifical 7-point stencil.
-    if (i_r > 0) {
-        const int bottom_left = grid.index(i_r - 1, i_theta_M1);
-        const int top_left    = grid.index(i_r - 1, i_theta_P1);
-        value += -0.25 * (art[left] + art[bottom]) * x[bottom_left]; /* Bottom Left */
-        value += +0.25 * (art[left] + art[top]) * x[top_left]; /* Top Left */
+    if ((i_r == 0 && DirBC_Interior) || (i_r == grid.nr() - 1)) {
+        const int center = grid.index(i_r, i_theta);
+        result[center]   = x[center];
+        return;
     }
 
+    const int i_theta_M1     = grid.wrapThetaIndex(i_theta - 1);
+    const int i_theta_P1     = grid.wrapThetaIndex(i_theta + 1);
+    const int i_theta_Across = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2);
+
+    const int left         = grid.index(i_r, i_theta_Across);
+    const int bottom       = grid.index(i_r, i_theta_M1);
+    const int center       = grid.index(i_r, i_theta);
+    const int top          = grid.index(i_r, i_theta_P1);
     const int bottom_right = grid.index(i_r + 1, i_theta_M1);
+    const int right        = grid.index(i_r + 1, i_theta);
     const int top_right    = grid.index(i_r + 1, i_theta_P1);
-    value += +0.25 * (art[right] + art[bottom]) * x[bottom_right]; /* Bottom Right */
-    value += -0.25 * (art[right] + art[top]) * x[top_right]; /* Top Right */
 
-    result[center] = value;
+    // Across origin: h1 gets replaced with 2 * R0.
+    const double h1 = 2.0 * grid.radius(0);
+    const double h2 = grid.radialSpacing(i_r);
+    const double k1 = grid.angularSpacing(i_theta - 1);
+    const double k2 = grid.angularSpacing(i_theta);
+
+    const double coeff1 = 0.5 * (k1 + k2) / h1;
+    const double coeff2 = 0.5 * (k1 + k2) / h2;
+    const double coeff3 = 0.5 * (h1 + h2) / k1;
+    const double coeff4 = 0.5 * (h1 + h2) / k2;
+    const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2);
+
+    result[center] = (+coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) * x[center] /* beta_{i,j} */
+
+                      + coeff1 * (arr[center] + arr[left]) * (x[center] - x[left]) /* Center: (Left) - Left */
+                      + coeff2 * (arr[center] + arr[right]) * (x[center] - x[right]) /* Center: (Right) - Right */
+                      + coeff3 * (att[center] + att[bottom]) * (x[center] - x[bottom]) /* Center: (Bottom) - Bottom */
+                      + coeff4 * (att[center] + att[top]) * (x[center] - x[top]) /* Center: (Top) - Top */
+
+                      + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
+                      - (art[right] + art[top]) * x[top_right] /* Top Right */
+    );
 }
 
 } // namespace residual_take
@@ -67,7 +104,8 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_take(const int i_r, const int i_
 template <class LevelCacheType>
 void ResidualTake<LevelCacheType>::applySystemOperator(Vector<double> result, ConstVector<double> x) const
 {
-    using residual_take::node_apply_a_take;
+    using residual_take::applySystemMatrixTakeBoundary;
+    using residual_take::applySystemMatrixTakeInterior;
 
     assert(result.size() == x.size());
 
@@ -88,26 +126,50 @@ void ResidualTake<LevelCacheType>::applySystemOperator(Vector<double> result, Co
 
     // The For loop matches circular access pattern */
     Kokkos::parallel_for(
-        "Residual Take: Apply System Operator (Circular)",
+        "Residual Take: Apply System Operator Boundary (Circular)",
         Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
             {0, 0}, // Starting point of the index space
+            {1, grid.ntheta()} // Ending point of the index space
+            ),
+        // Kokkos lambda function to execute for each point in the index space
+        KOKKOS_LAMBDA(const int i_r, const int i_theta) {
+            applySystemMatrixTakeBoundary(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF,
+                                          coeff_beta);
+        });
+
+    Kokkos::parallel_for(
+        "Residual Take: Apply System Operator Interior (Circular)",
+        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+            {1, 0}, // Starting point of the index space
             {grid.numberSmootherCircles(), grid.ntheta()} // Ending point of the index space
             ),
         // Kokkos lambda function to execute for each point in the index space
         KOKKOS_LAMBDA(const int i_r, const int i_theta) {
-            node_apply_a_take(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, coeff_beta);
+            applySystemMatrixTakeInterior(i_r, i_theta, grid, result, x, arr, att, art, detDF, coeff_beta);
         });
 
     /* For loop matches radial access pattern */
     Kokkos::parallel_for(
-        "Residual Take: Apply System Operator (Radial)",
+        "Residual Take: Apply System Operator Interior (Radial)",
         Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
             {0, grid.numberSmootherCircles()}, // Starting point of the index space
+            {grid.ntheta(), grid.nr() - 1} // Ending point of the index space
+            ),
+        // Kokkos lambda function to execute for each point in the index space
+        KOKKOS_LAMBDA(const int i_theta, const int i_r) {
+            applySystemMatrixTakeInterior(i_r, i_theta, grid, result, x, arr, att, art, detDF, coeff_beta);
+        });
+
+    Kokkos::parallel_for(
+        "Residual Take: Apply System Operator Boundary (Radial)",
+        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+            {0, grid.nr() - 1}, // Starting point of the index space
             {grid.ntheta(), grid.nr()} // Ending point of the index space
             ),
         // Kokkos lambda function to execute for each point in the index space
         KOKKOS_LAMBDA(const int i_theta, const int i_r) {
-            node_apply_a_take(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, coeff_beta);
+            applySystemMatrixTakeBoundary(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF,
+                                          coeff_beta);
         });
 
     Kokkos::fence();
diff --git a/include/Residual/ResidualTake/residualTake.h b/include/Residual/ResidualTake/residualTake.h
index c6577590..baad5398 100644
--- a/include/Residual/ResidualTake/residualTake.h
+++ b/include/Residual/ResidualTake/residualTake.h
@@ -10,8 +10,6 @@ class ResidualTake : public Residual<LevelCacheType>
 {
 public:
     explicit ResidualTake(const PolarGrid& grid, const LevelCacheType& level_cache, const bool DirBC_Interior);
-    KOKKOS_DEFAULTED_FUNCTION ResidualTake(const ResidualTake&) = default;
-    KOKKOS_DEFAULTED_FUNCTION ~ResidualTake() override          = default;
 
     void applySystemOperator(Vector<double> result, ConstVector<double> x) const final;
     void computeResidual(Vector<double> result, ConstVector<double> rhs, ConstVector<double> x) const final;
diff --git a/include/Smoother/SmootherGive/applyAscOrtho.inl b/include/Smoother/SmootherGive/applyAscOrtho.inl
index 9e0b29a5..fda04684 100644
--- a/include/Smoother/SmootherGive/applyAscOrtho.inl
+++ b/include/Smoother/SmootherGive/applyAscOrtho.inl
@@ -40,16 +40,16 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveInside(int i_r, in
         const double coeff1 = 0.5 * (k1 + k2) / h1;
         const int left      = grid.index(i_r - 1, i_theta);
         result[center] -= (-coeff1 * arr * x[left]); /* Left */
-        result[bottom] -= (+0.25 * art * x[left]); /* Top Left */
-        result[top] -= (-0.25 * art * x[left]); /* Bottom Left */
+        result[bottom] -= (+art * x[left]); /* Top Left */
+        result[top] -= (-art * x[left]); /* Bottom Left */
     }
 
     const double h2     = grid.radialSpacing(i_r);
     const double coeff2 = 0.5 * (k1 + k2) / h2;
     const int right     = grid.index(i_r + 1, i_theta);
     result[center] -= (-coeff2 * arr * x[right]); /* Right */
-    result[bottom] -= (-0.25 * art * x[right]); /* Top Right */
-    result[top] -= (+0.25 * art * x[right]); /* Bottom Right */
+    result[bottom] -= (-art * x[right]); /* Top Right */
+    result[top] -= (+art * x[right]); /* Bottom Right */
 }
 
 template <class LevelCacheType>
@@ -85,8 +85,8 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveOutside(int i_r, i
         const double coeff1 = 0.5 * (k1 + k2) / h1;
         const int left      = grid.index(i_r - 1, i_theta);
         result[left] -= (-coeff1 * arr * x[center] /* Right */
-                         - 0.25 * art * x[top] /* Top Right */
-                         + 0.25 * art * x[bottom]); /* Bottom Right */
+                         - art * x[top] /* Top Right */
+                         + art * x[bottom]); /* Bottom Right */
     }
     /* Fill result(i+1,j) */
     if (i_r < grid.numberSmootherCircles() - 1) {
@@ -94,8 +94,8 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveOutside(int i_r, i
         const double coeff2 = 0.5 * (k1 + k2) / h2;
         const int right     = grid.index(i_r + 1, i_theta);
         result[right] -= (-coeff2 * arr * x[center] /* Left */
-                          + 0.25 * art * x[top] /* Top Left */
-                          - 0.25 * art * x[bottom]); /* Bottom Left */
+                          + art * x[top] /* Top Left */
+                          - art * x[bottom]); /* Bottom Left */
     }
 }
 
@@ -138,14 +138,14 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialGiveInside(int i_r, in
 
     if (grid.numberSmootherCircles() <= i_r && i_r <= grid.nr() - 1) {
         const int left = grid.index(i_r - 1, i_theta);
-        result[left] -= (-0.25 * art * x[top] /* Top Right */
-                         + 0.25 * art * x[bottom]); /* Bottom Right */
+        result[left] -= (-art * x[top] /* Top Right */
+                         + art * x[bottom]); /* Bottom Right */
     }
 
     if (grid.numberSmootherCircles() - 1 <= i_r && i_r < grid.nr() - 2) {
         const int right = grid.index(i_r + 1, i_theta);
-        result[right] -= (+0.25 * art * x[top] /* Top Left    */
-                          - 0.25 * art * x[bottom]); /* Bottom Left */
+        result[right] -= (+art * x[top] /* Top Left    */
+                          - art * x[bottom]); /* Bottom Left */
     }
 
     if (i_r == grid.numberSmootherCircles() - 1) {
@@ -214,12 +214,12 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialGiveOutside(int i_r, i
 
         /* Fill result(i,j-1) */
         result[bottom] -= (-coeff3 * att * x[center] /* Top */
-                           - 0.25 * art * x[right] /* Top Right */
-                           + 0.25 * art * x[left]); /* Top Left */
+                           - art * x[right] /* Top Right */
+                           + art * x[left]); /* Top Left */
         /* Fill result(i,j+1) */
         result[top] -= (-coeff4 * att * x[center] /* Bottom */
-                        + 0.25 * art * x[right] /* Bottom Right */
-                        - 0.25 * art * x[left]); /* Bottom Left */
+                        + art * x[right] /* Bottom Right */
+                        - art * x[left]); /* Bottom Left */
     }
 }
 
diff --git a/include/Smoother/SmootherTake/applyAscOrtho.inl b/include/Smoother/SmootherTake/applyAscOrtho.inl
index d8b0af86..df7fa8ea 100644
--- a/include/Smoother/SmootherTake/applyAscOrtho.inl
+++ b/include/Smoother/SmootherTake/applyAscOrtho.inl
@@ -3,101 +3,159 @@
 namespace smoother_take
 {
 
-static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleTake(const int i_r, const int i_theta, const PolarGrid& grid,
-                                                               bool DirBC_Interior, ConstVector<double>& x,
-                                                               ConstVector<double>& rhs, Vector<double>& result,
-                                                               ConstVector<double>& arr, ConstVector<double>& att,
-                                                               ConstVector<double>& art)
+static KOKKOS_INLINE_FUNCTION void
+nodeApplyAscOrthoCircleTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, ConstVector<double>& x,
+                                    ConstVector<double>& rhs, Vector<double>& result, ConstVector<double>& arr,
+                                    ConstVector<double>& att, ConstVector<double>& art)
 {
-    KOKKOS_ASSERT(i_r >= 0 && i_r < grid.numberSmootherCircles());
+    KOKKOS_ASSERT(0 < i_r && i_r < grid.numberSmootherCircles());
 
-    const int center = grid.index(i_r, i_theta);
-
-    result[center] = rhs[center];
+    const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
+    const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
 
-    if (i_r == 0 && DirBC_Interior)
-        return;
+    const int bottom_left  = grid.index(i_r - 1, i_theta_M1);
+    const int left         = grid.index(i_r - 1, i_theta);
+    const int top_left     = grid.index(i_r - 1, i_theta_P1);
+    const int bottom       = grid.index(i_r, i_theta_M1);
+    const int center       = grid.index(i_r, i_theta);
+    const int top          = grid.index(i_r, i_theta_P1);
+    const int bottom_right = grid.index(i_r + 1, i_theta_M1);
+    const int right        = grid.index(i_r + 1, i_theta);
+    const int top_right    = grid.index(i_r + 1, i_theta_P1);
 
-    const double k1 = grid.angularSpacing(i_theta - 1);
+    const double h1 = grid.radialSpacing(i_r - 1);
+    const double h2 = grid.radialSpacing(i_r);
+    const double k1 = grid.angularSpacing(i_theta_M1);
     const double 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);
+    const double coeff1 = 0.5 * (k1 + k2) / h1;
+    const double coeff2 = 0.5 * (k1 + k2) / h2;
 
-    const int bottom = grid.index(i_r, i_theta_M1);
-    const int top    = grid.index(i_r, i_theta_P1);
+    result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */
+                                    - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */
 
-    if (i_r > 0) {
-        const double h1     = grid.radialSpacing(i_r - 1);
-        const double coeff1 = 0.5 * (k1 + k2) / h1;
+                                    - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
+                                    + (art[left] + art[top]) * x[top_left] /* Top Left */
+                                    + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
+                                    - (art[right] + art[top]) * x[top_right] /* Top Right */
+                                   );
+}
 
-        const int bottom_left = grid.index(i_r - 1, i_theta_M1);
-        const int left        = grid.index(i_r - 1, i_theta);
-        const int top_left    = grid.index(i_r - 1, i_theta_P1);
+static KOKKOS_INLINE_FUNCTION void
+nodeApplyAscOrthoCircleTakeBoundary(const int i_r, const int i_theta, const PolarGrid& grid, bool DirBC_Interior,
+                                    ConstVector<double>& x, ConstVector<double>& rhs, Vector<double>& result,
+                                    ConstVector<double>& arr, ConstVector<double>& att, ConstVector<double>& art)
+{
+    KOKKOS_ASSERT(i_r == 0);
 
-        result[center] -= (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */
-                           - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
-                           + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */
-        );
+    if (DirBC_Interior) {
+        const int center = grid.index(i_r, i_theta);
+        result[center]   = rhs[center];
+        return;
     }
 
-    const double h2     = grid.radialSpacing(i_r);
-    const double coeff2 = 0.5 * (k1 + k2) / h2;
+    const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
+    const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
 
+    const int bottom       = grid.index(i_r, i_theta_M1);
+    const int center       = grid.index(i_r, i_theta);
+    const int top          = grid.index(i_r, i_theta_P1);
     const int bottom_right = grid.index(i_r + 1, i_theta_M1);
     const int right        = grid.index(i_r + 1, i_theta);
     const int top_right    = grid.index(i_r + 1, i_theta_P1);
 
-    result[center] -= (-coeff2 * (arr[center] + arr[right]) * x[right] /* Right */
-                       + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
-                       - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */
+    const double h2 = grid.radialSpacing(i_r);
+    const double k1 = grid.angularSpacing(i_theta_M1);
+    const double k2 = grid.angularSpacing(i_theta);
 
-    );
+    const double coeff2 = 0.5 * (k1 + k2) / h2;
+
+    result[center] = rhs[center] - (-coeff2 * (arr[center] + arr[right]) * x[right] /* Right */
+                                    + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
+                                    - (art[right] + art[top]) * x[top_right] /* Top Right */
+                                   );
 }
 
-static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(const int i_r, const int i_theta, const PolarGrid& grid,
-                                                               bool DirBC_Interior, ConstVector<double>& x,
-                                                               ConstVector<double>& rhs, Vector<double>& result,
-                                                               ConstVector<double>& arr, ConstVector<double>& att,
-                                                               ConstVector<double>& art)
+static KOKKOS_INLINE_FUNCTION void
+nodeApplyAscOrthoRadialTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, ConstVector<double>& x,
+                                    ConstVector<double>& rhs, Vector<double>& result, ConstVector<double>& arr,
+                                    ConstVector<double>& att, ConstVector<double>& art)
 {
-    KOKKOS_ASSERT(i_r >= grid.numberSmootherCircles() && i_r < grid.nr());
-
-    const int center = grid.index(i_r, i_theta);
+    KOKKOS_ASSERT(grid.numberSmootherCircles() < i_r && i_r < grid.nr() - 2);
 
-    result[center] = rhs[center];
+    const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
+    const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
 
-    if (i_r == grid.nr() - 1)
-        return;
+    const int bottom_left  = grid.index(i_r - 1, i_theta_M1);
+    const int bottom       = grid.index(i_r, i_theta_M1);
+    const int bottom_right = grid.index(i_r + 1, i_theta_M1);
+    const int left         = grid.index(i_r - 1, i_theta);
+    const int center       = grid.index(i_r, i_theta);
+    const int right        = grid.index(i_r + 1, i_theta);
+    const int top_left     = grid.index(i_r - 1, i_theta_P1);
+    const int top          = grid.index(i_r, i_theta_P1);
+    const int top_right    = grid.index(i_r + 1, i_theta_P1);
 
     const double h1 = grid.radialSpacing(i_r - 1);
     const double h2 = grid.radialSpacing(i_r);
-    const double k1 = grid.angularSpacing(i_theta - 1);
+    const double k1 = grid.angularSpacing(i_theta_M1);
     const double k2 = grid.angularSpacing(i_theta);
 
     const double coeff3 = 0.5 * (h1 + h2) / k1;
     const double coeff4 = 0.5 * (h1 + h2) / k2;
 
+    result[center] = rhs[center] - (-coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
+                                    - coeff4 * (att[center] + att[top]) * x[top] /* Top */
+
+                                    - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
+                                    + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
+                                    + (art[left] + art[top]) * x[top_left] /* Top Left */
+                                    - (art[right] + art[top]) * x[top_right] /* Top Right */
+                                   );
+}
+
+static KOKKOS_INLINE_FUNCTION void
+nodeApplyAscOrthoRadialTakeBoundary(const int i_r, const int i_theta, const PolarGrid& grid, ConstVector<double>& x,
+                                    ConstVector<double>& rhs, Vector<double>& result, ConstVector<double>& arr,
+                                    ConstVector<double>& att, ConstVector<double>& art)
+{
+    KOKKOS_ASSERT(i_r == grid.numberSmootherCircles() || i_r == grid.nr() - 2 || i_r == grid.nr() - 1);
+
+    if (i_r == grid.nr() - 1) {
+        const int center = grid.index(i_r, i_theta);
+        result[center]   = rhs[center];
+        return;
+    }
+
     const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
     const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
 
     const int bottom_left  = grid.index(i_r - 1, i_theta_M1);
-    const int left         = grid.index(i_r - 1, i_theta);
-    const int top_left     = grid.index(i_r - 1, i_theta_P1);
     const int bottom       = grid.index(i_r, i_theta_M1);
-    const int top          = grid.index(i_r, i_theta_P1);
     const int bottom_right = grid.index(i_r + 1, i_theta_M1);
+    const int left         = grid.index(i_r - 1, i_theta);
+    const int center       = grid.index(i_r, i_theta);
     const int right        = grid.index(i_r + 1, i_theta);
+    const int top_left     = grid.index(i_r - 1, i_theta_P1);
+    const int top          = grid.index(i_r, i_theta_P1);
     const int top_right    = grid.index(i_r + 1, i_theta_P1);
 
-    result[center] -= (-coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
-                       - coeff4 * (att[center] + att[top]) * x[top] /* Top */
+    const double h1 = grid.radialSpacing(i_r - 1);
+    const double h2 = grid.radialSpacing(i_r);
+    const double k1 = grid.angularSpacing(i_theta_M1);
+    const double k2 = grid.angularSpacing(i_theta);
 
-                       - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
-                       + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
-                       + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */
-                       - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */
-    );
+    const double coeff3 = 0.5 * (h1 + h2) / k1;
+    const double coeff4 = 0.5 * (h1 + h2) / k2;
+
+    result[center] = rhs[center] - (-coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */
+                                    - coeff4 * (att[center] + att[top]) * x[top] /* Top */
+
+                                    - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
+                                    + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
+                                    + (art[left] + art[top]) * x[top_left] /* Top Left */
+                                    - (art[right] + art[top]) * x[top_right] /* Top Right */
+                                   );
 
     if (i_r == grid.numberSmootherCircles()) {
         const double coeff1 = 0.5 * (k1 + k2) / h1;
@@ -120,7 +178,8 @@ template <class LevelCacheType>
 void SmootherTake<LevelCacheType>::applyAscOrthoBlackCircleSection(ConstVector<double> x, ConstVector<double> rhs,
                                                                    Vector<double> temp)
 {
-    using smoother_take::nodeApplyAscOrthoCircleTake;
+    using smoother_take::nodeApplyAscOrthoCircleTakeBoundary;
+    using smoother_take::nodeApplyAscOrthoCircleTakeInterior;
 
     const PolarGrid& grid             = Smoother<LevelCacheType>::grid_;
     const LevelCacheType& level_cache = Smoother<LevelCacheType>::level_cache_;
@@ -136,17 +195,44 @@ void SmootherTake<LevelCacheType>::applyAscOrthoBlackCircleSection(ConstVector<d
     const int start_black_circles = (grid.numberSmootherCircles() % 2 == 0) ? 1 : 0;
     const int num_black_circles   = (grid.numberSmootherCircles() - start_black_circles + 1) / 2;
 
-    Kokkos::parallel_for(
-        "Smoother Take: ApplyAscOrtho (Black Circular)",
-        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
-            {0, 0}, // Starting point of the index space
-            {num_black_circles, grid.ntheta()} // Ending point of the index space
-            ),
-        // Kokkos lambda function to execute for each point in the index space
-        KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
-            int i_r = start_black_circles + circle_task * 2;
-            nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art);
-        });
+    if (start_black_circles == 0) {
+        Kokkos::parallel_for(
+            "Smoother Take: ApplyAscOrtho (Black Circular)",
+            Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+                {0, 0}, // Starting point of the index space
+                {1, grid.ntheta()} // Ending point of the index space
+                ),
+            // Kokkos lambda function to execute for each point in the index space
+            KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
+                int i_r = start_black_circles + circle_task * 2;
+                nodeApplyAscOrthoCircleTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art);
+            });
+
+        Kokkos::parallel_for(
+            "Smoother Take: ApplyAscOrtho (Black Circular)",
+            Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+                {1, 0}, // Starting point of the index space
+                {num_black_circles, grid.ntheta()} // Ending point of the index space
+                ),
+            // Kokkos lambda function to execute for each point in the index space
+            KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
+                int i_r = start_black_circles + circle_task * 2;
+                nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+            });
+    }
+    else {
+        Kokkos::parallel_for(
+            "Smoother Take: ApplyAscOrtho (Black Circular)",
+            Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+                {0, 0}, // Starting point of the index space
+                {num_black_circles, grid.ntheta()} // Ending point of the index space
+                ),
+            // Kokkos lambda function to execute for each point in the index space
+            KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
+                int i_r = start_black_circles + circle_task * 2;
+                nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+            });
+    }
 
     Kokkos::fence();
 }
@@ -155,7 +241,8 @@ template <class LevelCacheType>
 void SmootherTake<LevelCacheType>::applyAscOrthoWhiteCircleSection(ConstVector<double> x, ConstVector<double> rhs,
                                                                    Vector<double> temp)
 {
-    using smoother_take::nodeApplyAscOrthoCircleTake;
+    using smoother_take::nodeApplyAscOrthoCircleTakeBoundary;
+    using smoother_take::nodeApplyAscOrthoCircleTakeInterior;
 
     const PolarGrid& grid             = Smoother<LevelCacheType>::grid_;
     const LevelCacheType& level_cache = Smoother<LevelCacheType>::level_cache_;
@@ -171,17 +258,44 @@ void SmootherTake<LevelCacheType>::applyAscOrthoWhiteCircleSection(ConstVector<d
     const int start_white_circles = (grid.numberSmootherCircles() % 2 == 0) ? 0 : 1;
     const int num_white_circles   = (grid.numberSmootherCircles() - start_white_circles + 1) / 2;
 
-    Kokkos::parallel_for(
-        "Smoother Take: ApplyAscOrtho (White Circular)",
-        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
-            {0, 0}, // Starting point of the index space
-            {num_white_circles, grid.ntheta()} // Ending point of the index space
-            ),
-        // Kokkos lambda function to execute for each point in the index space
-        KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
-            const int i_r = start_white_circles + circle_task * 2;
-            nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art);
-        });
+    if (start_white_circles == 0) {
+        Kokkos::parallel_for(
+            "Smoother Take: ApplyAscOrtho (White Circular)",
+            Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+                {0, 0}, // Starting point of the index space
+                {1, grid.ntheta()} // Ending point of the index space
+                ),
+            // Kokkos lambda function to execute for each point in the index space
+            KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
+                const int i_r = start_white_circles + circle_task * 2;
+                nodeApplyAscOrthoCircleTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art);
+            });
+
+        Kokkos::parallel_for(
+            "Smoother Take: ApplyAscOrtho (White Circular)",
+            Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+                {1, 0}, // Starting point of the index space
+                {num_white_circles, grid.ntheta()} // Ending point of the index space
+                ),
+            // Kokkos lambda function to execute for each point in the index space
+            KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
+                const int i_r = start_white_circles + circle_task * 2;
+                nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+            });
+    }
+    else {
+        Kokkos::parallel_for(
+            "Smoother Take: ApplyAscOrtho (White Circular)",
+            Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+                {0, 0}, // Starting point of the index space
+                {num_white_circles, grid.ntheta()} // Ending point of the index space
+                ),
+            // Kokkos lambda function to execute for each point in the index space
+            KOKKOS_LAMBDA(const int circle_task, const int i_theta) {
+                const int i_r = start_white_circles + circle_task * 2;
+                nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+            });
+    }
 
     Kokkos::fence();
 }
@@ -190,7 +304,8 @@ template <class LevelCacheType>
 void SmootherTake<LevelCacheType>::applyAscOrthoBlackRadialSection(ConstVector<double> x, ConstVector<double> rhs,
                                                                    Vector<double> temp)
 {
-    using smoother_take::nodeApplyAscOrthoRadialTake;
+    using smoother_take::nodeApplyAscOrthoRadialTakeBoundary;
+    using smoother_take::nodeApplyAscOrthoRadialTakeInterior;
 
     const PolarGrid& grid             = Smoother<LevelCacheType>::grid_;
     const LevelCacheType& level_cache = Smoother<LevelCacheType>::level_cache_;
@@ -210,12 +325,36 @@ void SmootherTake<LevelCacheType>::applyAscOrthoBlackRadialSection(ConstVector<d
         "Smoother Take: ApplyAscOrtho (Black Radial)",
         Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
             {0, grid.numberSmootherCircles()}, // Starting point of the index space
+            {num_black_radial_lines, grid.numberSmootherCircles() + 1} // Ending point of the index space
+            ),
+        // Kokkos lambda function to execute for each point in the index space
+        KOKKOS_LAMBDA(const int radial_task, const int i_r) {
+            const int i_theta = start_black_radials + radial_task * 2;
+            nodeApplyAscOrthoRadialTakeBoundary(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+        });
+
+    Kokkos::parallel_for(
+        "Smoother Take: ApplyAscOrtho (Black Radial)",
+        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+            {0, grid.numberSmootherCircles() + 1}, // Starting point of the index space
+            {num_black_radial_lines, grid.nr() - 2} // Ending point of the index space
+            ),
+        // Kokkos lambda function to execute for each point in the index space
+        KOKKOS_LAMBDA(const int radial_task, const int i_r) {
+            const int i_theta = start_black_radials + radial_task * 2;
+            nodeApplyAscOrthoRadialTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+        });
+
+    Kokkos::parallel_for(
+        "Smoother Take: ApplyAscOrtho (Black Radial)",
+        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+            {0, grid.nr() - 2}, // Starting point of the index space
             {num_black_radial_lines, grid.nr()} // Ending point of the index space
             ),
         // Kokkos lambda function to execute for each point in the index space
         KOKKOS_LAMBDA(const int radial_task, const int i_r) {
             const int i_theta = start_black_radials + radial_task * 2;
-            nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art);
+            nodeApplyAscOrthoRadialTakeBoundary(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
         });
 
     Kokkos::fence();
@@ -225,7 +364,8 @@ template <class LevelCacheType>
 void SmootherTake<LevelCacheType>::applyAscOrthoWhiteRadialSection(ConstVector<double> x, ConstVector<double> rhs,
                                                                    Vector<double> temp)
 {
-    using smoother_take::nodeApplyAscOrthoRadialTake;
+    using smoother_take::nodeApplyAscOrthoRadialTakeBoundary;
+    using smoother_take::nodeApplyAscOrthoRadialTakeInterior;
 
     const PolarGrid& grid             = Smoother<LevelCacheType>::grid_;
     const LevelCacheType& level_cache = Smoother<LevelCacheType>::level_cache_;
@@ -245,12 +385,36 @@ void SmootherTake<LevelCacheType>::applyAscOrthoWhiteRadialSection(ConstVector<d
         "Smoother Take: ApplyAscOrtho (White Radial)",
         Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
             {0, grid.numberSmootherCircles()}, // Starting point of the index space
+            {num_white_radial_lines, grid.numberSmootherCircles() + 1} // Ending point of the index space
+            ),
+        // Kokkos lambda function to execute for each point in the index space
+        KOKKOS_LAMBDA(const int radial_task, const int i_r) {
+            const int i_theta = start_white_radials + radial_task * 2;
+            nodeApplyAscOrthoRadialTakeBoundary(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+        });
+
+    Kokkos::parallel_for(
+        "Smoother Take: ApplyAscOrtho (White Radial)",
+        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+            {0, grid.numberSmootherCircles() + 1}, // Starting point of the index space
+            {num_white_radial_lines, grid.nr() - 2} // Ending point of the index space
+            ),
+        // Kokkos lambda function to execute for each point in the index space
+        KOKKOS_LAMBDA(const int radial_task, const int i_r) {
+            const int i_theta = start_white_radials + radial_task * 2;
+            nodeApplyAscOrthoRadialTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
+        });
+
+    Kokkos::parallel_for(
+        "Smoother Take: ApplyAscOrtho (White Radial)",
+        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>( // Rank of the index space
+            {0, grid.nr() - 2}, // Starting point of the index space
             {num_white_radial_lines, grid.nr()} // Ending point of the index space
             ),
         // Kokkos lambda function to execute for each point in the index space
         KOKKOS_LAMBDA(const int radial_task, const int i_r) {
             const int i_theta = start_white_radials + radial_task * 2;
-            nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art);
+            nodeApplyAscOrthoRadialTakeBoundary(i_r, i_theta, grid, x, rhs, temp, arr, att, art);
         });
 
     Kokkos::fence();
diff --git a/include/Smoother/SmootherTake/smootherTake.h b/include/Smoother/SmootherTake/smootherTake.h
index d04ee584..2e137899 100644
--- a/include/Smoother/SmootherTake/smootherTake.h
+++ b/include/Smoother/SmootherTake/smootherTake.h
@@ -55,8 +55,6 @@ class SmootherTake : public Smoother<LevelCacheType>
     // Builds the A_sc smoother matrices and prepares the solvers.
     explicit SmootherTake(const PolarGrid& grid, const LevelCacheType& level_cache, bool DirBC_Interior);
 
-    KOKKOS_DEFAULTED_FUNCTION SmootherTake(const SmootherTake&) = default;
-
     // Performs one full coupled smoothing sweep:
     //   BC -> WC -> BR -> WR
     // using temp as RHS workspace.
diff --git a/scripts/tutorial/run.sh b/scripts/tutorial/run.sh
index fb9ea8a5..572d4dd4 100755
--- a/scripts/tutorial/run.sh
+++ b/scripts/tutorial/run.sh
@@ -168,4 +168,4 @@ export OMP_NUM_THREADS=$maxOpenMPThreads
     --maxIterations $maxIterations \
     --residualNormType $residualNormType \
     --absoluteTolerance $absoluteTolerance \
-    --relativeTolerance $relativeTolerance
+    --relativeTolerance $relativeTolerance
\ No newline at end of file