diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h
index 2df76754..7b4fccff 100644
--- a/include/GMGPolar/gmgpolar.h
+++ b/include/GMGPolar/gmgpolar.h
@@ -176,7 +176,7 @@ class GMGPolar : public IGMGPolar
     void updateResidualNorms(Level<DomainGeometry, DensityProfileCoefficients>& level, int iteration,
                              double& initial_residual_norm, double& current_residual_norm,
                              double& current_relative_residual_norm);
-    void initRhsHierarchy(Vector<double> rhs);
+    void initRhsHierarchy();
     void applyMultigridIterations(Level<DomainGeometry, DensityProfileCoefficients>& level, MultigridCycleType cycle,
                                   int iterations);
     void applyExtrapolatedMultigridIterations(Level<DomainGeometry, DensityProfileCoefficients>& level,
diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h
index 559f25e4..7fd3601c 100644
--- a/include/GMGPolar/solver.h
+++ b/include/GMGPolar/solver.h
@@ -318,7 +318,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::solvePCG(double& init
 
     // z = M^{-1} * r (preconditioned residual)
     if (PCG_FMG_) {
-        initRhsHierarchy(level.rhs());
+        initRhsHierarchy();
         fullMultigridApproximation(PCG_FMG_cycle_, PCG_FMG_iterations_);
     }
     else {
@@ -392,7 +392,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::solvePCG(double& init
 
         // z = M^{-1} * r (preconditioned residual)
         if (PCG_FMG_) {
-            initRhsHierarchy(level.rhs());
+            initRhsHierarchy();
             fullMultigridApproximation(PCG_FMG_cycle_, PCG_FMG_iterations_);
         }
         else {
@@ -577,9 +577,8 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::applyExtrapolation(in
 // =============================================================================
 
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
-void GMGPolar<DomainGeometry, DensityProfileCoefficients>::initRhsHierarchy(Vector<double> rhs)
+void GMGPolar<DomainGeometry, DensityProfileCoefficients>::initRhsHierarchy()
 {
-    Kokkos::deep_copy(levels_[0].rhs(), rhs);
     for (int level_depth = 0; level_depth < number_of_levels_ - 1; ++level_depth) {
         Level<DomainGeometry, DensityProfileCoefficients>& current_level = levels_[level_depth];
         Level<DomainGeometry, DensityProfileCoefficients>& next_level    = levels_[level_depth + 1];