diff --git a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h
index 403dc95b..2e249315 100644
--- a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h
+++ b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h
@@ -1,7 +1,7 @@
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
 void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigrid_F_Cycle(int level_depth,
                                                                                           Vector<double> solution,
-                                                                                          Vector<double> rhs,
+                                                                                          ConstVector<double> rhs,
                                                                                           Vector<double> residual)
 {
     assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
@@ -94,8 +94,14 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
         assign(next_level.residual(), 0.0);
 
         /* Step 3: Solve for the error by recursively calling the multigrid cycle. */
-        multigrid_F_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
-        multigrid_V_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
+        multigrid_F_Cycle(level_depth + 1,
+                          next_level.residual(), // error (solution)
+                          next_level.error_correction(), // coarse residual (rhs)
+                          next_level.solution()); // workspace
+        multigrid_V_Cycle(level_depth + 1,
+                          next_level.residual(), // error (solution)
+                          next_level.error_correction(), // coarse residual (rhs)
+                          next_level.solution()); // workspace
     }
 
     /* Interpolate the correction */
diff --git a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h
index 2462f398..1877c621 100644
--- a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h
+++ b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h
@@ -1,7 +1,7 @@
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
 void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigrid_V_Cycle(int level_depth,
                                                                                           Vector<double> solution,
-                                                                                          Vector<double> rhs,
+                                                                                          ConstVector<double> rhs,
                                                                                           Vector<double> residual)
 {
     assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
@@ -94,7 +94,10 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
         assign(next_level.residual(), 0.0);
 
         /* Step 3: Solve for the error by recursively calling the multigrid cycle. */
-        multigrid_V_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
+        multigrid_V_Cycle(level_depth + 1,
+                          next_level.residual(), // error (solution)
+                          next_level.error_correction(), // coarse residual (rhs)
+                          next_level.solution()); // workspace
     }
 
     /* Interpolate the correction */
diff --git a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h
index 55dc257e..173ed826 100644
--- a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h
+++ b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h
@@ -1,7 +1,7 @@
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
 void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigrid_W_Cycle(int level_depth,
                                                                                           Vector<double> solution,
-                                                                                          Vector<double> rhs,
+                                                                                          ConstVector<double> rhs,
                                                                                           Vector<double> residual)
 {
     assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
@@ -94,8 +94,14 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
         assign(next_level.residual(), 0.0);
 
         /* Step 3: Solve for the error by recursively calling the multigrid cycle. */
-        multigrid_W_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
-        multigrid_W_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
+        multigrid_W_Cycle(level_depth + 1,
+                          next_level.residual(), // error (solution)
+                          next_level.error_correction(), // coarse residual (rhs)
+                          next_level.solution()); // workspace
+        multigrid_W_Cycle(level_depth + 1,
+                          next_level.residual(), // error (solution)
+                          next_level.error_correction(), // coarse residual (rhs)
+                          next_level.solution()); // workspace
     }
 
     /* Interpolate the correction */
diff --git a/include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h b/include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h
index 913bd972..a2eaa57b 100644
--- a/include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h
+++ b/include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h
@@ -1,63 +1,71 @@
+#pragma once
+
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
 void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_F_Cycle(int level_depth, Vector<double> solution,
-                                                                             Vector<double> rhs,
+                                                                             ConstVector<double> rhs,
                                                                              Vector<double> residual)
 {
-    assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
+    assert(0 <= level_depth && level_depth < number_of_levels_);
 
     std::chrono::high_resolution_clock::time_point start_MGC;
     if (level_depth == 0) {
         start_MGC = std::chrono::high_resolution_clock::now();
     }
 
-    Level<DomainGeometry>& level      = levels_[level_depth];
-    Level<DomainGeometry>& next_level = levels_[level_depth + 1];
-
-    auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
-
-    /* ------------ */
-    /* Presmoothing */
-    for (int i = 0; i < pre_smoothing_steps_; i++) {
-        level.smoothing(solution, rhs, residual);
-    }
-
-    auto end_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
-    t_avg_MGC_preSmoothing_ += std::chrono::duration<double>(end_MGC_preSmoothing - start_MGC_preSmoothing).count();
-
-    /* ---------------------- */
-    /* Coarse grid correction */
-    /* ---------------------- */
-
-    auto start_MGC_residual = std::chrono::high_resolution_clock::now();
-
-    /* Compute the residual */
-    level.computeResidual(residual, rhs, solution);
+    /* ------------------------ */
+    /* Solve A * solution = rhs */
+    /* ------------------------ */
+    if (level_depth == number_of_levels_ - 1) {
+        /* ------------------------------ */
+        /* Coarsest level: solve directly */
+        /* ------------------------------ */
+        Level<DomainGeometry>& level = levels_[level_depth];
 
-    auto end_MGC_residual = std::chrono::high_resolution_clock::now();
-    t_avg_MGC_residual_ += std::chrono::duration<double>(end_MGC_residual - start_MGC_residual).count();
+        /* Step 1: Copy rhs in solution */
+        Kokkos::deep_copy(solution, rhs);
 
-    /* -------------------------- */
-    /* Solve A * error = residual */
-    if (level_depth + 1 == number_of_levels_ - 1) {
-        /* --------------------- */
-        /* Using a direct solver */
-        /* --------------------- */
-
-        /* Step 1: Restrict the residual */
-        restriction(level_depth, next_level.residual(), residual);
-
-        /* Step 2: Solve for the error in place */
+        /* Step 2: Solve for the solution in place */
         auto start_MGC_directSolver = std::chrono::high_resolution_clock::now();
 
-        next_level.directSolveInPlace(next_level.residual());
+        level.directSolveInPlace(solution);
 
         auto end_MGC_directSolver = std::chrono::high_resolution_clock::now();
         t_avg_MGC_directSolver_ += std::chrono::duration<double>(end_MGC_directSolver - start_MGC_directSolver).count();
     }
     else {
-        /* ------------------------------------------ */
+        /* ------------------------------------------------------- */
+        /* Multigrid V-cycle on current level:                     */
+        /* presmoothing, coarse-grid correction, and postsmoothing */
+        /* ------------------------------------------------------- */
+        Level<DomainGeometry>& level      = levels_[level_depth];
+        Level<DomainGeometry>& next_level = levels_[level_depth + 1];
+
+        auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
+
+        /* ------------ */
+        /* Presmoothing */
+        for (int i = 0; i < pre_smoothing_steps_; i++) {
+            level.smoothing(solution, rhs, residual);
+        }
+
+        auto end_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
+        t_avg_MGC_preSmoothing_ += std::chrono::duration<double>(end_MGC_preSmoothing - start_MGC_preSmoothing).count();
+
+        /* ----------------------------- */
+        /* Compute fine-grid residual    */
+        /* residual = rhs - A * solution */
+        /* ----------------------------- */
+        auto start_MGC_residual = std::chrono::high_resolution_clock::now();
+
+        level.computeResidual(residual, rhs, solution);
+
+        auto end_MGC_residual = std::chrono::high_resolution_clock::now();
+        t_avg_MGC_residual_ += std::chrono::duration<double>(end_MGC_residual - start_MGC_residual).count();
+
+        /* -------------------------- */
+        /* Solve A * error = residual */
+        /* -------------------------- */
         /* By recursively calling the multigrid cycle */
-        /* ------------------------------------------ */
 
         /* Step 1: Restrict the residual. */
         restriction(level_depth, next_level.error_correction(), residual);
@@ -66,29 +74,36 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_F_Cycle(int
         assign(next_level.residual(), 0.0);
 
         /* Step 3: Solve for the error by recursively calling the multigrid cycle. */
-        multigrid_F_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
-        multigrid_V_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
-    }
-
-    /* Interpolate the correction */
-    prolongation(level_depth + 1, residual, next_level.residual());
-
-    /* Compute the corrected approximation: u = u + error */
-    add(solution, ConstVector<double>(residual));
-
-    auto start_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
-
-    /* ------------- */
-    /* Postsmoothing */
-    for (int i = 0; i < post_smoothing_steps_; i++) {
-        level.smoothing(solution, rhs, residual);
+        multigrid_F_Cycle(level_depth + 1,
+                          next_level.residual(), // error (solution)
+                          next_level.error_correction(), // coarse residual (rhs)
+                          next_level.solution()); // workspace
+        multigrid_V_Cycle(level_depth + 1,
+                          next_level.residual(), // error (solution)
+                          next_level.error_correction(), // coarse residual (rhs)
+                          next_level.solution()); // workspace
+
+        /* Interpolate the correction */
+        prolongation(level_depth + 1, residual, next_level.residual());
+
+        /* Compute the corrected approximation: u = u + error */
+        add(solution, ConstVector<double>(residual));
+
+        auto start_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
+
+        /* ------------- */
+        /* Postsmoothing */
+        for (int i = 0; i < post_smoothing_steps_; i++) {
+            level.smoothing(solution, rhs, residual);
+        }
+
+        auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
+        t_avg_MGC_postSmoothing_ +=
+            std::chrono::duration<double>(end_MGC_postSmoothing - start_MGC_postSmoothing).count();
     }
 
-    auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
-    t_avg_MGC_postSmoothing_ += std::chrono::duration<double>(end_MGC_postSmoothing - start_MGC_postSmoothing).count();
-
     if (level_depth == 0) {
         auto end_MGC = std::chrono::high_resolution_clock::now();
         t_avg_MGC_total_ += std::chrono::duration<double>(end_MGC - start_MGC).count();
     }
-}
+}
\ No newline at end of file
diff --git a/include/GMGPolar/MultigridMethods/multigrid_V_Cycle.h b/include/GMGPolar/MultigridMethods/multigrid_V_Cycle.h
index 931ec0c7..596928b3 100644
--- a/include/GMGPolar/MultigridMethods/multigrid_V_Cycle.h
+++ b/include/GMGPolar/MultigridMethods/multigrid_V_Cycle.h
@@ -1,63 +1,71 @@
+#pragma once
+
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
 void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_V_Cycle(int level_depth, Vector<double> solution,
-                                                                             Vector<double> rhs,
+                                                                             ConstVector<double> rhs,
                                                                              Vector<double> residual)
 {
-    assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
+    assert(0 <= level_depth && level_depth < number_of_levels_);
 
     std::chrono::high_resolution_clock::time_point start_MGC;
     if (level_depth == 0) {
         start_MGC = std::chrono::high_resolution_clock::now();
     }
 
-    Level<DomainGeometry>& level      = levels_[level_depth];
-    Level<DomainGeometry>& next_level = levels_[level_depth + 1];
-
-    auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
-
-    /* ------------ */
-    /* Presmoothing */
-    for (int i = 0; i < pre_smoothing_steps_; i++) {
-        level.smoothing(solution, rhs, residual);
-    }
-
-    auto end_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
-    t_avg_MGC_preSmoothing_ += std::chrono::duration<double>(end_MGC_preSmoothing - start_MGC_preSmoothing).count();
-
-    /* ---------------------- */
-    /* Coarse grid correction */
-    /* ---------------------- */
-
-    auto start_MGC_residual = std::chrono::high_resolution_clock::now();
+    /* ------------------------ */
+    /* Solve A * solution = rhs */
+    /* ------------------------ */
+    if (level_depth == number_of_levels_ - 1) {
+        /* ------------------------------ */
+        /* Coarsest level: solve directly */
+        /* ------------------------------ */
+        Level<DomainGeometry>& level = levels_[level_depth];
 
-    /* Compute the residual */
-    level.computeResidual(residual, rhs, solution);
+        /* Step 1: Copy rhs in solution */
+        Kokkos::deep_copy(solution, rhs);
 
-    auto end_MGC_residual = std::chrono::high_resolution_clock::now();
-    t_avg_MGC_residual_ += std::chrono::duration<double>(end_MGC_residual - start_MGC_residual).count();
-
-    /* -------------------------- */
-    /* Solve A * error = residual */
-    if (level_depth + 1 == number_of_levels_ - 1) {
-        /* --------------------- */
-        /* Using a direct solver */
-        /* --------------------- */
-
-        /* Step 1: Restrict the residual */
-        restriction(level_depth, next_level.residual(), residual);
-
-        /* Step 2: Solve for the error in place */
+        /* Step 2: Solve for the solution in place */
         auto start_MGC_directSolver = std::chrono::high_resolution_clock::now();
 
-        next_level.directSolveInPlace(next_level.residual());
+        level.directSolveInPlace(solution);
 
         auto end_MGC_directSolver = std::chrono::high_resolution_clock::now();
         t_avg_MGC_directSolver_ += std::chrono::duration<double>(end_MGC_directSolver - start_MGC_directSolver).count();
     }
     else {
-        /* ------------------------------------------ */
+        /* ------------------------------------------------------- */
+        /* Multigrid V-cycle on current level:                     */
+        /* presmoothing, coarse-grid correction, and postsmoothing */
+        /* ------------------------------------------------------- */
+        Level<DomainGeometry>& level      = levels_[level_depth];
+        Level<DomainGeometry>& next_level = levels_[level_depth + 1];
+
+        auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
+
+        /* ------------ */
+        /* Presmoothing */
+        for (int i = 0; i < pre_smoothing_steps_; i++) {
+            level.smoothing(solution, rhs, residual);
+        }
+
+        auto end_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
+        t_avg_MGC_preSmoothing_ += std::chrono::duration<double>(end_MGC_preSmoothing - start_MGC_preSmoothing).count();
+
+        /* ----------------------------- */
+        /* Compute fine-grid residual    */
+        /* residual = rhs - A * solution */
+        /* ----------------------------- */
+        auto start_MGC_residual = std::chrono::high_resolution_clock::now();
+
+        level.computeResidual(residual, rhs, solution);
+
+        auto end_MGC_residual = std::chrono::high_resolution_clock::now();
+        t_avg_MGC_residual_ += std::chrono::duration<double>(end_MGC_residual - start_MGC_residual).count();
+
+        /* -------------------------- */
+        /* Solve A * error = residual */
+        /* -------------------------- */
         /* By recursively calling the multigrid cycle */
-        /* ------------------------------------------ */
 
         /* Step 1: Restrict the residual. */
         restriction(level_depth, next_level.error_correction(), residual);
@@ -66,28 +74,32 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_V_Cycle(int
         assign(next_level.residual(), 0.0);
 
         /* Step 3: Solve for the error by recursively calling the multigrid cycle. */
-        multigrid_V_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
-    }
+        multigrid_V_Cycle(level_depth + 1,
+                          next_level.residual(), // error (solution)
+                          next_level.error_correction(), // coarse residual (rhs)
+                          next_level.solution()); // workspace
 
-    /* Interpolate the correction */
-    prolongation(level_depth + 1, residual, next_level.residual());
+        /* Interpolate the correction */
+        prolongation(level_depth + 1, residual, next_level.residual());
 
-    /* Compute the corrected approximation: u = u + error */
-    add(solution, ConstVector<double>(residual));
+        /* Compute the corrected approximation: u = u + error */
+        add(solution, ConstVector<double>(residual));
 
-    auto start_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
+        auto start_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
 
-    /* ------------- */
-    /* Postsmoothing */
-    for (int i = 0; i < post_smoothing_steps_; i++) {
-        level.smoothing(solution, rhs, residual);
-    }
+        /* ------------- */
+        /* Postsmoothing */
+        for (int i = 0; i < post_smoothing_steps_; i++) {
+            level.smoothing(solution, rhs, residual);
+        }
 
-    auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
-    t_avg_MGC_postSmoothing_ += std::chrono::duration<double>(end_MGC_postSmoothing - start_MGC_postSmoothing).count();
+        auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
+        t_avg_MGC_postSmoothing_ +=
+            std::chrono::duration<double>(end_MGC_postSmoothing - start_MGC_postSmoothing).count();
+    }
 
     if (level_depth == 0) {
         auto end_MGC = std::chrono::high_resolution_clock::now();
         t_avg_MGC_total_ += std::chrono::duration<double>(end_MGC - start_MGC).count();
     }
-}
+}
\ No newline at end of file
diff --git a/include/GMGPolar/MultigridMethods/multigrid_W_Cycle.h b/include/GMGPolar/MultigridMethods/multigrid_W_Cycle.h
index a2d1c522..1c52d56d 100644
--- a/include/GMGPolar/MultigridMethods/multigrid_W_Cycle.h
+++ b/include/GMGPolar/MultigridMethods/multigrid_W_Cycle.h
@@ -1,63 +1,71 @@
+#pragma once
+
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
 void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_W_Cycle(int level_depth, Vector<double> solution,
-                                                                             Vector<double> rhs,
+                                                                             ConstVector<double> rhs,
                                                                              Vector<double> residual)
 {
-    assert(0 <= level_depth && level_depth < number_of_levels_ - 1);
+    assert(0 <= level_depth && level_depth < number_of_levels_);
 
     std::chrono::high_resolution_clock::time_point start_MGC;
     if (level_depth == 0) {
         start_MGC = std::chrono::high_resolution_clock::now();
     }
 
-    Level<DomainGeometry>& level      = levels_[level_depth];
-    Level<DomainGeometry>& next_level = levels_[level_depth + 1];
-
-    auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
-
-    /* ------------ */
-    /* Presmoothing */
-    for (int i = 0; i < pre_smoothing_steps_; i++) {
-        level.smoothing(solution, rhs, residual);
-    }
-
-    auto end_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
-    t_avg_MGC_preSmoothing_ += std::chrono::duration<double>(end_MGC_preSmoothing - start_MGC_preSmoothing).count();
-
-    /* ---------------------- */
-    /* Coarse grid correction */
-    /* ---------------------- */
-
-    auto start_MGC_residual = std::chrono::high_resolution_clock::now();
-
-    /* Compute the residual */
-    level.computeResidual(residual, rhs, solution);
+    /* ------------------------ */
+    /* Solve A * solution = rhs */
+    /* ------------------------ */
+    if (level_depth == number_of_levels_ - 1) {
+        /* ------------------------------ */
+        /* Coarsest level: solve directly */
+        /* ------------------------------ */
+        Level<DomainGeometry>& level = levels_[level_depth];
 
-    auto end_MGC_residual = std::chrono::high_resolution_clock::now();
-    t_avg_MGC_residual_ += std::chrono::duration<double>(end_MGC_residual - start_MGC_residual).count();
+        /* Step 1: Copy rhs in solution */
+        Kokkos::deep_copy(solution, rhs);
 
-    /* -------------------------- */
-    /* Solve A * error = residual */
-    if (level_depth + 1 == number_of_levels_ - 1) {
-        /* --------------------- */
-        /* Using a direct solver */
-        /* --------------------- */
-
-        /* Step 1: Restrict the residual */
-        restriction(level_depth, next_level.residual(), residual);
-
-        /* Step 2: Solve for the error in place */
+        /* Step 2: Solve for the solution in place */
         auto start_MGC_directSolver = std::chrono::high_resolution_clock::now();
 
-        next_level.directSolveInPlace(next_level.residual());
+        level.directSolveInPlace(solution);
 
         auto end_MGC_directSolver = std::chrono::high_resolution_clock::now();
         t_avg_MGC_directSolver_ += std::chrono::duration<double>(end_MGC_directSolver - start_MGC_directSolver).count();
     }
     else {
-        /* ------------------------------------------ */
+        /* ------------------------------------------------------- */
+        /* Multigrid V-cycle on current level:                     */
+        /* presmoothing, coarse-grid correction, and postsmoothing */
+        /* ------------------------------------------------------- */
+        Level<DomainGeometry>& level      = levels_[level_depth];
+        Level<DomainGeometry>& next_level = levels_[level_depth + 1];
+
+        auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
+
+        /* ------------ */
+        /* Presmoothing */
+        for (int i = 0; i < pre_smoothing_steps_; i++) {
+            level.smoothing(solution, rhs, residual);
+        }
+
+        auto end_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
+        t_avg_MGC_preSmoothing_ += std::chrono::duration<double>(end_MGC_preSmoothing - start_MGC_preSmoothing).count();
+
+        /* ----------------------------- */
+        /* Compute fine-grid residual    */
+        /* residual = rhs - A * solution */
+        /* ----------------------------- */
+        auto start_MGC_residual = std::chrono::high_resolution_clock::now();
+
+        level.computeResidual(residual, rhs, solution);
+
+        auto end_MGC_residual = std::chrono::high_resolution_clock::now();
+        t_avg_MGC_residual_ += std::chrono::duration<double>(end_MGC_residual - start_MGC_residual).count();
+
+        /* -------------------------- */
+        /* Solve A * error = residual */
+        /* -------------------------- */
         /* By recursively calling the multigrid cycle */
-        /* ------------------------------------------ */
 
         /* Step 1: Restrict the residual. */
         restriction(level_depth, next_level.error_correction(), residual);
@@ -66,29 +74,36 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_W_Cycle(int
         assign(next_level.residual(), 0.0);
 
         /* Step 3: Solve for the error by recursively calling the multigrid cycle. */
-        multigrid_W_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
-        multigrid_W_Cycle(level_depth + 1, next_level.residual(), next_level.error_correction(), next_level.solution());
-    }
-
-    /* Interpolate the correction */
-    prolongation(level_depth + 1, residual, next_level.residual());
-
-    /* Compute the corrected approximation: u = u + error */
-    add(solution, ConstVector<double>(residual));
-
-    auto start_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
-
-    /* ------------- */
-    /* Postsmoothing */
-    for (int i = 0; i < post_smoothing_steps_; i++) {
-        level.smoothing(solution, rhs, residual);
+        multigrid_W_Cycle(level_depth + 1,
+                          next_level.residual(), // error (solution)
+                          next_level.error_correction(), // coarse residual (rhs)
+                          next_level.solution()); // workspace
+        multigrid_W_Cycle(level_depth + 1,
+                          next_level.residual(), // error (solution)
+                          next_level.error_correction(), // coarse residual (rhs)
+                          next_level.solution()); // workspace
+
+        /* Interpolate the correction */
+        prolongation(level_depth + 1, residual, next_level.residual());
+
+        /* Compute the corrected approximation: u = u + error */
+        add(solution, ConstVector<double>(residual));
+
+        auto start_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
+
+        /* ------------- */
+        /* Postsmoothing */
+        for (int i = 0; i < post_smoothing_steps_; i++) {
+            level.smoothing(solution, rhs, residual);
+        }
+
+        auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
+        t_avg_MGC_postSmoothing_ +=
+            std::chrono::duration<double>(end_MGC_postSmoothing - start_MGC_postSmoothing).count();
     }
 
-    auto end_MGC_postSmoothing = std::chrono::high_resolution_clock::now();
-    t_avg_MGC_postSmoothing_ += std::chrono::duration<double>(end_MGC_postSmoothing - start_MGC_postSmoothing).count();
-
     if (level_depth == 0) {
         auto end_MGC = std::chrono::high_resolution_clock::now();
         t_avg_MGC_total_ += std::chrono::duration<double>(end_MGC - start_MGC).count();
     }
-}
+}
\ No newline at end of file
diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h
index cbd070ae..4d5ee2c6 100644
--- a/include/GMGPolar/gmgpolar.h
+++ b/include/GMGPolar/gmgpolar.h
@@ -167,14 +167,14 @@ class GMGPolar : public IGMGPolar
 
     /* ------------------- */
     /* Multigrid Functions */
-    void multigrid_V_Cycle(int level_depth, Vector<double> solution, Vector<double> rhs, Vector<double> residual);
-    void multigrid_W_Cycle(int level_depth, Vector<double> solution, Vector<double> rhs, Vector<double> residual);
-    void multigrid_F_Cycle(int level_depth, Vector<double> solution, Vector<double> rhs, Vector<double> residual);
-    void extrapolated_multigrid_V_Cycle(int level_depth, Vector<double> solution, Vector<double> rhs,
+    void multigrid_V_Cycle(int level_depth, Vector<double> solution, ConstVector<double> rhs, Vector<double> residual);
+    void multigrid_W_Cycle(int level_depth, Vector<double> solution, ConstVector<double> rhs, Vector<double> residual);
+    void multigrid_F_Cycle(int level_depth, Vector<double> solution, ConstVector<double> rhs, Vector<double> residual);
+    void extrapolated_multigrid_V_Cycle(int level_depth, Vector<double> solution, ConstVector<double> rhs,
                                         Vector<double> residual);
-    void extrapolated_multigrid_W_Cycle(int level_depth, Vector<double> solution, Vector<double> rhs,
+    void extrapolated_multigrid_W_Cycle(int level_depth, Vector<double> solution, ConstVector<double> rhs,
                                         Vector<double> residual);
-    void extrapolated_multigrid_F_Cycle(int level_depth, Vector<double> solution, Vector<double> rhs,
+    void extrapolated_multigrid_F_Cycle(int level_depth, Vector<double> solution, ConstVector<double> rhs,
                                         Vector<double> residual);
 
     /* ----------------------- */
diff --git a/include/GMGPolar/setup.h b/include/GMGPolar/setup.h
index d9194b21..24f6ca29 100644
--- a/include/GMGPolar/setup.h
+++ b/include/GMGPolar/setup.h
@@ -174,7 +174,7 @@ int GMGPolar<DomainGeometry, DensityProfileCoefficients>::chooseNumberOfLevels(c
     const int minAngularDivisions = 4;
 
     // Minimum level for Multigrid
-    const int multigridMinLevel = 2;
+    const int multigridMinLevel = (extrapolation_ == ExtrapolationType::NONE) ? 1 : 2;
 
     // Calculate radial maximum level
     int radialNodes    = finestGrid.nr();
diff --git a/tests/GMGPolar/solve_tests.cpp b/tests/GMGPolar/solve_tests.cpp
index e6f69c1c..aa6c53e3 100644
--- a/tests/GMGPolar/solve_tests.cpp
+++ b/tests/GMGPolar/solve_tests.cpp
@@ -545,6 +545,108 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<double, 8e-5>, // expected_l2_error
         std::integral_constant<double, 3e-4>, // expected_inf_error
         std::integral_constant<double, 0.7> // expected_residual_reduction
+    >,
+    std::tuple<
+        CzarnyGeometry,
+        SonnendruckerGyroCoefficients,
+        PolarR6_Boundary_CzarnyGeometry,
+        PolarR6_Sonnendrucker_CzarnyGeometry,
+        PolarR6_CzarnyGeometry,
+        std::integral_constant<double, 1e-8>, // R0
+        std::integral_constant<double, 1.2>, // Rmax
+        std::integral_constant<int, 4>, // nrExp
+        std::integral_constant<int, 5>, // nthetaExp
+        std::integral_constant<double, 0.66 * 1.2>, // refinementRadius
+        std::integral_constant<int, 0>, // anisotropicFactor
+        std::integral_constant<int, 0>, // divideBy2
+        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // maxOpenMPThreads
+        std::integral_constant<bool, false>, // DirBC_Interior
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
+        std::integral_constant<bool, false>, // cacheDomainGeometry
+        std::integral_constant<ExtrapolationType, ExtrapolationType::NONE>, // extrapolation
+        std::integral_constant<int, 1>, // maxLevels
+        std::integral_constant<MultigridCycleType, MultigridCycleType::V_CYCLE>, // multigridCycle
+        std::integral_constant<bool, false>, // FMG
+        std::integral_constant<int, 1>, // FMG_iterations
+        std::integral_constant<MultigridCycleType, MultigridCycleType::F_CYCLE>, // FMG_cycle
+        std::integral_constant<int, 50>, // maxIterations
+        std::integral_constant<ResidualNormType, ResidualNormType::EUCLIDEAN>, // residualNormType
+        std::integral_constant<double, 1e-12>, // absoluteTolerance
+        std::integral_constant<double, 1e-12>, // relativeTolerance
+        std::integral_constant<int, 1>, // expected_iterations
+        std::integral_constant<double, 2e-5>, // expected_l2_error
+        std::integral_constant<double, 5e-5>, // expected_inf_error
+        std::integral_constant<double, 0.7> // expected_residual_reduction
+    >,
+    std::tuple<
+        CzarnyGeometry,
+        SonnendruckerGyroCoefficients,
+        PolarR6_Boundary_CzarnyGeometry,
+        PolarR6_Sonnendrucker_CzarnyGeometry,
+        PolarR6_CzarnyGeometry,
+        std::integral_constant<double, 1e-8>, // R0
+        std::integral_constant<double, 1.2>, // Rmax
+        std::integral_constant<int, 4>, // nrExp
+        std::integral_constant<int, 5>, // nthetaExp
+        std::integral_constant<double, 0.66 * 1.2>, // refinementRadius
+        std::integral_constant<int, 0>, // anisotropicFactor
+        std::integral_constant<int, 0>, // divideBy2
+        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // maxOpenMPThreads
+        std::integral_constant<bool, false>, // DirBC_Interior
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
+        std::integral_constant<bool, false>, // cacheDomainGeometry
+        std::integral_constant<ExtrapolationType, ExtrapolationType::NONE>, // extrapolation
+        std::integral_constant<int, 1>, // maxLevels
+        std::integral_constant<MultigridCycleType, MultigridCycleType::F_CYCLE>, // multigridCycle
+        std::integral_constant<bool, false>, // FMG
+        std::integral_constant<int, 1>, // FMG_iterations
+        std::integral_constant<MultigridCycleType, MultigridCycleType::F_CYCLE>, // FMG_cycle
+        std::integral_constant<int, 50>, // maxIterations
+        std::integral_constant<ResidualNormType, ResidualNormType::EUCLIDEAN>, // residualNormType
+        std::integral_constant<double, 1e-12>, // absoluteTolerance
+        std::integral_constant<double, 1e-12>, // relativeTolerance
+        std::integral_constant<int, 1>, // expected_iterations
+        std::integral_constant<double, 2e-5>, // expected_l2_error
+        std::integral_constant<double, 5e-5>, // expected_inf_error
+        std::integral_constant<double, 0.7> // expected_residual_reduction
+    >,
+    std::tuple<
+        CzarnyGeometry,
+        SonnendruckerGyroCoefficients,
+        PolarR6_Boundary_CzarnyGeometry,
+        PolarR6_Sonnendrucker_CzarnyGeometry,
+        PolarR6_CzarnyGeometry,
+        std::integral_constant<double, 1e-8>, // R0
+        std::integral_constant<double, 1.2>, // Rmax
+        std::integral_constant<int, 4>, // nrExp
+        std::integral_constant<int, 5>, // nthetaExp
+        std::integral_constant<double, 0.66 * 1.2>, // refinementRadius
+        std::integral_constant<int, 0>, // anisotropicFactor
+        std::integral_constant<int, 0>, // divideBy2
+        std::integral_constant<int, 0>, // verbose
+        std::integral_constant<int, 1>, // maxOpenMPThreads
+        std::integral_constant<bool, false>, // DirBC_Interior
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
+        std::integral_constant<bool, false>, // cacheDomainGeometry
+        std::integral_constant<ExtrapolationType, ExtrapolationType::NONE>, // extrapolation
+        std::integral_constant<int, 1>, // maxLevels
+        std::integral_constant<MultigridCycleType, MultigridCycleType::W_CYCLE>, // multigridCycle
+        std::integral_constant<bool, false>, // FMG
+        std::integral_constant<int, 1>, // FMG_iterations
+        std::integral_constant<MultigridCycleType, MultigridCycleType::F_CYCLE>, // FMG_cycle
+        std::integral_constant<int, 50>, // maxIterations
+        std::integral_constant<ResidualNormType, ResidualNormType::EUCLIDEAN>, // residualNormType
+        std::integral_constant<double, 1e-12>, // absoluteTolerance
+        std::integral_constant<double, 1e-12>, // relativeTolerance
+        std::integral_constant<int, 1>, // expected_iterations
+        std::integral_constant<double, 2e-5>, // expected_l2_error
+        std::integral_constant<double, 5e-5>, // expected_inf_error
+        std::integral_constant<double, 0.7> // expected_residual_reduction
     >
 >;
 // clang-format on