diff --git a/include/Definitions/global_definitions.h b/include/Definitions/global_definitions.h
index d90ae9a6d..379b32166 100644
--- a/include/Definitions/global_definitions.h
+++ b/include/Definitions/global_definitions.h
@@ -6,8 +6,8 @@
 
 enum class StencilDistributionMethod
 {
-    CPU_TAKE = 0,
-    CPU_GIVE = 1
+    TAKE = 0,
+    GIVE = 1
 };
 
 /* Multigrid Cycle Types */
@@ -34,13 +34,6 @@ enum class ExtrapolationType
     COMBINED                     = 3,
 };
 
-/* Smoother Colors */
-enum class SmootherColor
-{
-    Black = 0,
-    White = 1,
-};
-
 /* -----------*/
 /* Test Cases */
 /* -----------*/
diff --git a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h
index b36e98a41..f20524608 100644
--- a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h
+++ b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h
@@ -17,8 +17,8 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
     /* ------------------------------ */
     /* Extrapolated multigrid F-cycle */
     /* ------------------------------ */
-    Level<DomainGeometry, DensityProfileCoefficients>& level      = levels_[level_depth];
-    Level<DomainGeometry, DensityProfileCoefficients>& next_level = levels_[level_depth + 1];
+    LevelType& level      = levels_[level_depth];
+    LevelType& next_level = levels_[level_depth + 1];
 
     /* ------------ */
     /* Presmoothing */
diff --git a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h
index 36e6c67eb..d7d12b165 100644
--- a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h
+++ b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h
@@ -17,8 +17,8 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
     /* ------------------------------ */
     /* Extrapolated multigrid V-cycle */
     /* ------------------------------ */
-    Level<DomainGeometry, DensityProfileCoefficients>& level      = levels_[level_depth];
-    Level<DomainGeometry, DensityProfileCoefficients>& next_level = levels_[level_depth + 1];
+    LevelType& level      = levels_[level_depth];
+    LevelType& next_level = levels_[level_depth + 1];
 
     /* ------------ */
     /* Presmoothing */
diff --git a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h
index 0b1577352..9a07efa3b 100644
--- a/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h
+++ b/include/GMGPolar/MultigridMethods/extrapolated_multigrid_W_Cycle.h
@@ -17,8 +17,8 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
     /* ------------------------------ */
     /* Extrapolated multigrid W-cycle */
     /* ------------------------------ */
-    Level<DomainGeometry, DensityProfileCoefficients>& level      = levels_[level_depth];
-    Level<DomainGeometry, DensityProfileCoefficients>& next_level = levels_[level_depth + 1];
+    LevelType& level      = levels_[level_depth];
+    LevelType& next_level = levels_[level_depth + 1];
 
     /* ------------ */
     /* Presmoothing */
diff --git a/include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h b/include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h
index abeede937..304a6de9d 100644
--- a/include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h
+++ b/include/GMGPolar/MultigridMethods/multigrid_F_Cycle.h
@@ -16,7 +16,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_F_Cycle(int
         /* ---------------------------------------------------- */
         /* Coarsest level: solve A * x = rhs using DirectSolver */
         /* ---------------------------------------------------- */
-        Level<DomainGeometry, DensityProfileCoefficients>& coarsest_level = levels_[level_depth];
+        LevelType& coarsest_level = levels_[level_depth];
 
         /* Step 1: Copy rhs in solution */
         Kokkos::deep_copy(solution, rhs);
@@ -33,8 +33,8 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_F_Cycle(int
         /* ----------------- */
         /* Multigrid F-cycle */
         /* ----------------- */
-        Level<DomainGeometry, DensityProfileCoefficients>& level      = levels_[level_depth];
-        Level<DomainGeometry, DensityProfileCoefficients>& next_level = levels_[level_depth + 1];
+        LevelType& level      = levels_[level_depth];
+        LevelType& next_level = levels_[level_depth + 1];
 
         /* ------------ */
         /* Presmoothing */
diff --git a/include/GMGPolar/MultigridMethods/multigrid_V_Cycle.h b/include/GMGPolar/MultigridMethods/multigrid_V_Cycle.h
index f8d6e1287..69c2ccbdf 100644
--- a/include/GMGPolar/MultigridMethods/multigrid_V_Cycle.h
+++ b/include/GMGPolar/MultigridMethods/multigrid_V_Cycle.h
@@ -16,7 +16,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_V_Cycle(int
         /* ---------------------------------------------------- */
         /* Coarsest level: solve A * x = rhs using DirectSolver */
         /* ---------------------------------------------------- */
-        Level<DomainGeometry, DensityProfileCoefficients>& coarsest_level = levels_[level_depth];
+        LevelType& coarsest_level = levels_[level_depth];
 
         /* Step 1: Copy rhs in solution */
         Kokkos::deep_copy(solution, rhs);
@@ -33,8 +33,8 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_V_Cycle(int
         /* ----------------- */
         /* Multigrid V-cycle */
         /* ----------------- */
-        Level<DomainGeometry, DensityProfileCoefficients>& level      = levels_[level_depth];
-        Level<DomainGeometry, DensityProfileCoefficients>& next_level = levels_[level_depth + 1];
+        LevelType& level      = levels_[level_depth];
+        LevelType& next_level = levels_[level_depth + 1];
 
         /* ------------ */
         /* Presmoothing */
diff --git a/include/GMGPolar/MultigridMethods/multigrid_W_Cycle.h b/include/GMGPolar/MultigridMethods/multigrid_W_Cycle.h
index cf59b1209..4dc233cfe 100644
--- a/include/GMGPolar/MultigridMethods/multigrid_W_Cycle.h
+++ b/include/GMGPolar/MultigridMethods/multigrid_W_Cycle.h
@@ -16,7 +16,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_W_Cycle(int
         /* ---------------------------------------------------- */
         /* Coarsest level: solve A * x = rhs using DirectSolver */
         /* ---------------------------------------------------- */
-        Level<DomainGeometry, DensityProfileCoefficients>& coarsest_level = levels_[level_depth];
+        LevelType& coarsest_level = levels_[level_depth];
 
         /* Step 1: Copy rhs in solution */
         Kokkos::deep_copy(solution, rhs);
@@ -33,8 +33,8 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::multigrid_W_Cycle(int
         /* ----------------- */
         /* Multigrid W-cycle */
         /* ----------------- */
-        Level<DomainGeometry, DensityProfileCoefficients>& level      = levels_[level_depth];
-        Level<DomainGeometry, DensityProfileCoefficients>& next_level = levels_[level_depth + 1];
+        LevelType& level      = levels_[level_depth];
+        LevelType& next_level = levels_[level_depth + 1];
 
         /* ------------ */
         /* Presmoothing */
diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h
index dc80f18ec..daed475ac 100644
--- a/include/GMGPolar/gmgpolar.h
+++ b/include/GMGPolar/gmgpolar.h
@@ -18,6 +18,8 @@ template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoeff
 class GMGPolar : public IGMGPolar
 {
 public:
+    using LevelType = Level<DomainGeometry, DensityProfileCoefficients>;
+
     /* ---------------------------------------------------------------------- */
     /* Constructor & Initialization                                           */
     /* ---------------------------------------------------------------------- */
@@ -94,7 +96,7 @@ class GMGPolar : public IGMGPolar
     /* ---------------- */
     /* Multigrid levels */
     int number_of_levels_;
-    std::vector<Level<DomainGeometry, DensityProfileCoefficients>> levels_;
+    std::vector<LevelType> levels_;
 
     /* ---------------------- */
     /* Interpolation operator */
@@ -145,9 +147,9 @@ class GMGPolar : public IGMGPolar
     /* --------------- */
     /* Setup Functions */
     template <concepts::BoundaryConditions BoundaryConditions, concepts::SourceTerm SourceTerm>
-    void build_rhs_f(const Level<DomainGeometry, DensityProfileCoefficients>& level, Vector<double> rhs_f,
-                     const BoundaryConditions& boundary_conditions, const SourceTerm& source_term);
-    void discretize_rhs_f(const Level<DomainGeometry, DensityProfileCoefficients>& level, Vector<double> rhs_f);
+    void build_rhs_f(const LevelType& level, Vector<double> rhs_f, const BoundaryConditions& boundary_conditions,
+                     const SourceTerm& source_term);
+    void discretize_rhs_f(const LevelType& level, Vector<double> rhs_f);
 
     /* --------------- */
     /* Solve Functions */
@@ -165,19 +167,13 @@ class GMGPolar : public IGMGPolar
     void solveMultigrid(double& initial_residual_norm, double& current_residual_norm,
                         double& current_relative_residual_norm);
     void solvePCG(double& initial_residual_norm, double& current_residual_norm, double& current_relative_residual_norm);
-    double residualNorm(const ResidualNormType& norm_type,
-                        const Level<DomainGeometry, DensityProfileCoefficients>& level,
-                        ConstVector<double> residual) const;
-    void evaluateExactError(Level<DomainGeometry, DensityProfileCoefficients>& level,
-                            HostConstVector<double> exact_solution);
-    void updateResidualNorms(Level<DomainGeometry, DensityProfileCoefficients>& level, int iteration,
-                             double& initial_residual_norm, double& current_residual_norm,
-                             double& current_relative_residual_norm);
+    double residualNorm(const ResidualNormType& norm_type, const LevelType& level, ConstVector<double> residual) const;
+    void evaluateExactError(LevelType& level, HostConstVector<double> exact_solution);
+    void updateResidualNorms(LevelType& level, int iteration, double& initial_residual_norm,
+                             double& current_residual_norm, double& current_relative_residual_norm);
     void initRhsHierarchy();
-    void applyMultigridIterations(Level<DomainGeometry, DensityProfileCoefficients>& level, MultigridCycleType cycle,
-                                  int iterations);
-    void applyExtrapolatedMultigridIterations(Level<DomainGeometry, DensityProfileCoefficients>& level,
-                                              MultigridCycleType cycle, int iterations);
+    void applyMultigridIterations(LevelType& level, MultigridCycleType cycle, int iterations);
+    void applyExtrapolatedMultigridIterations(LevelType& level, MultigridCycleType cycle, int iterations);
 
     /* ----------------- */
     /* Print information */
@@ -213,8 +209,7 @@ class GMGPolar : public IGMGPolar
     void writeToVTK(const std::filesystem::path& file_path, const PolarGrid& grid);
 
 private:
-    void writeToVTK(const std::filesystem::path& file_path,
-                    const Level<DomainGeometry, DensityProfileCoefficients>& level,
+    void writeToVTK(const std::filesystem::path& file_path, const LevelType& level,
                     HostConstVector<double> grid_function);
 };
 
diff --git a/include/GMGPolar/setup.h b/include/GMGPolar/setup.h
index f4e031b9a..abcdcb144 100644
--- a/include/GMGPolar/setup.h
+++ b/include/GMGPolar/setup.h
@@ -10,7 +10,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::setup()
 
     auto start_setup_createLevels = std::chrono::high_resolution_clock::now();
 
-    if (stencil_distribution_method_ == StencilDistributionMethod::CPU_TAKE) {
+    if (stencil_distribution_method_ == StencilDistributionMethod::TAKE) {
         if (!cache_density_profile_coefficients_ || !cache_domain_geometry_) {
             throw std::runtime_error("Error: Caching must be enabled for both density profile coefficients and domain "
                                      "geometry in 'Take' implementation strategy.");
@@ -190,8 +190,8 @@ int GMGPolar<DomainGeometry, DensityProfileCoefficients>::chooseNumberOfLevels(c
 }
 
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
-void GMGPolar<DomainGeometry, DensityProfileCoefficients>::discretize_rhs_f(
-    const Level<DomainGeometry, DensityProfileCoefficients>& level, Vector<double> rhs_f)
+void GMGPolar<DomainGeometry, DensityProfileCoefficients>::discretize_rhs_f(const LevelType& level,
+                                                                            Vector<double> rhs_f)
 {
     const PolarGrid& grid = level.grid();
     assert(std::ssize(rhs_f) == grid.numberOfNodes());
@@ -331,9 +331,9 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::discretize_rhs_f(
 
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
 template <concepts::BoundaryConditions BoundaryConditions, concepts::SourceTerm SourceTerm>
-void GMGPolar<DomainGeometry, DensityProfileCoefficients>::build_rhs_f(
-    const Level<DomainGeometry, DensityProfileCoefficients>& level, Vector<double> rhs_f,
-    const BoundaryConditions& boundary_conditions, const SourceTerm& source_term)
+void GMGPolar<DomainGeometry, DensityProfileCoefficients>::build_rhs_f(const LevelType& level, Vector<double> rhs_f,
+                                                                       const BoundaryConditions& boundary_conditions,
+                                                                       const SourceTerm& source_term)
 {
     const PolarGrid& grid(level.grid());
     assert(std::ssize(rhs_f) == grid.numberOfNodes());
@@ -424,7 +424,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::printSettings(const P
         std::cout << "Across the origin (Interior boundary condition)\n";
     }
 
-    if (stencil_distribution_method_ == StencilDistributionMethod::CPU_TAKE) {
+    if (stencil_distribution_method_ == StencilDistributionMethod::TAKE) {
         std::cout << "A-Take (Stencil Distribution)\n";
     }
     else {
diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h
index b7e9bdc28..63c708541 100644
--- a/include/GMGPolar/solver.h
+++ b/include/GMGPolar/solver.h
@@ -22,10 +22,10 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::solve(const BoundaryC
     int initial_rhs_f_levels = FMG_ ? number_of_levels_ : (extrapolation_ == ExtrapolationType::NONE ? 1 : 2);
     // Loop through the levels, injecting and discretizing rhs
     for (int level_depth = 0; level_depth < initial_rhs_f_levels; ++level_depth) {
-        Level<DomainGeometry, DensityProfileCoefficients>& current_level = levels_[level_depth];
+        LevelType& current_level = levels_[level_depth];
         // Inject rhs if there is a next level
         if (level_depth + 1 < initial_rhs_f_levels) {
-            Level<DomainGeometry, DensityProfileCoefficients>& next_level = levels_[level_depth + 1];
+            LevelType& next_level = levels_[level_depth + 1];
             injection(level_depth, next_level.rhs(), current_level.rhs());
         }
         // Discretize the rhs for the current level
@@ -82,7 +82,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::solve(const BoundaryC
     /* --------------------------------------- */
     /* Start Solver at finest level (depth 0)  */
     /* --------------------------------------- */
-    Level<DomainGeometry, DensityProfileCoefficients>& level = levels_[0];
+    LevelType& level = levels_[0];
 
     number_of_iterations_                 = 0;
     double initial_residual_norm          = 1.0;
@@ -201,8 +201,8 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::fullMultigridApproxim
                                                                                       int FMG_iterations)
 {
     // Start from the coarsest level
-    int coarsest_depth                                                = number_of_levels_ - 1;
-    Level<DomainGeometry, DensityProfileCoefficients>& coarsest_level = levels_[coarsest_depth];
+    int coarsest_depth        = number_of_levels_ - 1;
+    LevelType& coarsest_level = levels_[coarsest_depth];
 
     // Solve directly on the coarsest level
     Kokkos::deep_copy(coarsest_level.solution(), coarsest_level.rhs());
@@ -210,8 +210,8 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::fullMultigridApproxim
 
     // Prolongate the solution from the coarsest level up to the finest, while applying Multigrid Cycles on each level
     for (int depth = coarsest_depth; depth > 0; --depth) {
-        Level<DomainGeometry, DensityProfileCoefficients>& coarse_level = levels_[depth]; // Current coarse level
-        Level<DomainGeometry, DensityProfileCoefficients>& fine_level   = levels_[depth - 1]; // Next finer level
+        LevelType& coarse_level = levels_[depth]; // Current coarse level
+        LevelType& fine_level   = levels_[depth - 1]; // Next finer level
 
         // The bi-cubic FMG interpolation is of higher order
         FMGInterpolation(coarse_level.level_depth(), fine_level.solution(), coarse_level.solution());
@@ -233,7 +233,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::solveMultigrid(double
                                                                           double& current_residual_norm,
                                                                           double& current_relative_residual_norm)
 {
-    Level<DomainGeometry, DensityProfileCoefficients>& level = levels_[0];
+    LevelType& level = levels_[0];
 
     while (number_of_iterations_ < max_iterations_) {
         /* ----------------------- */
@@ -308,7 +308,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::solvePCG(double& init
                                                                     double& current_residual_norm,
                                                                     double& current_relative_residual_norm)
 {
-    Level<DomainGeometry, DensityProfileCoefficients>& level = levels_[0];
+    LevelType& level = levels_[0];
 
     // x = initial guess
     Kokkos::deep_copy(pcg_solution_, level.solution());
@@ -338,7 +338,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::solvePCG(double& init
         level.applySystemOperator(level.residual(), pcg_search_direction_);
         if (extrapolation_ != ExtrapolationType::NONE) {
             assert(number_of_levels_ > 1);
-            Level<DomainGeometry, DensityProfileCoefficients>& next_level = levels_[level.level_depth() + 1];
+            LevelType& next_level = levels_[level.level_depth() + 1];
             injection(0, next_level.solution(), pcg_search_direction_);
             next_level.applySystemOperator(next_level.residual(), next_level.solution());
             applyExtrapolation(0, level.residual(), next_level.residual());
@@ -421,8 +421,9 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::solvePCG(double& init
 // =============================================================================
 
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
-void GMGPolar<DomainGeometry, DensityProfileCoefficients>::applyMultigridIterations(
-    Level<DomainGeometry, DensityProfileCoefficients>& level, MultigridCycleType cycle, int iterations)
+void GMGPolar<DomainGeometry, DensityProfileCoefficients>::applyMultigridIterations(LevelType& level,
+                                                                                    MultigridCycleType cycle,
+                                                                                    int iterations)
 {
     for (int i = 0; i < iterations; i++) {
         switch (cycle) {
@@ -444,7 +445,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::applyMultigridIterati
 
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
 void GMGPolar<DomainGeometry, DensityProfileCoefficients>::applyExtrapolatedMultigridIterations(
-    Level<DomainGeometry, DensityProfileCoefficients>& level, MultigridCycleType cycle, int iterations)
+    LevelType& level, MultigridCycleType cycle, int iterations)
 {
     for (int i = 0; i < iterations; i++) {
         switch (cycle) {
@@ -469,13 +470,14 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::applyExtrapolatedMult
 // =============================================================================
 
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
-void GMGPolar<DomainGeometry, DensityProfileCoefficients>::updateResidualNorms(
-    Level<DomainGeometry, DensityProfileCoefficients>& level, int iteration, double& initial_residual_norm,
-    double& current_residual_norm, double& current_relative_residual_norm)
+void GMGPolar<DomainGeometry, DensityProfileCoefficients>::updateResidualNorms(LevelType& level, int iteration,
+                                                                               double& initial_residual_norm,
+                                                                               double& current_residual_norm,
+                                                                               double& current_relative_residual_norm)
 {
     level.computeResidual(level.residual(), level.rhs(), level.solution());
     if (extrapolation_ != ExtrapolationType::NONE) {
-        Level<DomainGeometry, DensityProfileCoefficients>& next_level = levels_[level.level_depth() + 1];
+        LevelType& next_level = levels_[level.level_depth() + 1];
         injection(level.level_depth(), next_level.solution(), level.solution());
         next_level.computeResidual(next_level.residual(), next_level.rhs(), next_level.solution());
         applyExtrapolation(level.level_depth(), level.residual(), next_level.residual());
@@ -504,9 +506,9 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::updateResidualNorms(
 }
 
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
-double GMGPolar<DomainGeometry, DensityProfileCoefficients>::residualNorm(
-    const ResidualNormType& norm_type, const Level<DomainGeometry, DensityProfileCoefficients>& level,
-    ConstVector<double> residual) const
+double GMGPolar<DomainGeometry, DensityProfileCoefficients>::residualNorm(const ResidualNormType& norm_type,
+                                                                          const LevelType& level,
+                                                                          ConstVector<double> residual) const
 {
     switch (norm_type) {
     case ResidualNormType::EUCLIDEAN:
@@ -579,8 +581,8 @@ template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoeff
 void GMGPolar<DomainGeometry, DensityProfileCoefficients>::initRhsHierarchy()
 {
     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];
+        LevelType& current_level = levels_[level_depth];
+        LevelType& next_level    = levels_[level_depth + 1];
         restriction(level_depth, next_level.rhs(), current_level.rhs());
     }
 }
diff --git a/include/GMGPolar/utils.h b/include/GMGPolar/utils.h
index df76f2973..1b0b5c60e 100644
--- a/include/GMGPolar/utils.h
+++ b/include/GMGPolar/utils.h
@@ -256,9 +256,9 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::writeToVTK(const std:
 }
 
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
-void GMGPolar<DomainGeometry, DensityProfileCoefficients>::writeToVTK(
-    const std::filesystem::path& file_path, const Level<DomainGeometry, DensityProfileCoefficients>& level,
-    HostConstVector<double> grid_function)
+void GMGPolar<DomainGeometry, DensityProfileCoefficients>::writeToVTK(const std::filesystem::path& file_path,
+                                                                      const LevelType& level,
+                                                                      HostConstVector<double> grid_function)
 {
     const PolarGrid& grid = level.grid();
 
diff --git a/include/Level/level.inl b/include/Level/level.inl
index 4b2832ac0..451487609 100644
--- a/include/Level/level.inl
+++ b/include/Level/level.inl
@@ -94,10 +94,10 @@ template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoeff
 void Level<DomainGeometry, DensityProfileCoefficients>::initializeResidual(
     const bool DirBC_Interior, const StencilDistributionMethod stencil_distribution_method)
 {
-    if (stencil_distribution_method == StencilDistributionMethod::CPU_TAKE) {
+    if (stencil_distribution_method == StencilDistributionMethod::TAKE) {
         op_residual_ = std::make_unique<ResidualTake<LevelCacheType>>(*grid_, *level_cache_, DirBC_Interior);
     }
-    else if (stencil_distribution_method == StencilDistributionMethod::CPU_GIVE) {
+    else if (stencil_distribution_method == StencilDistributionMethod::GIVE) {
         op_residual_ = std::make_unique<ResidualGive<LevelCacheType>>(*grid_, *level_cache_, DirBC_Interior);
     }
     if (!op_residual_)
@@ -127,10 +127,10 @@ template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoeff
 void Level<DomainGeometry, DensityProfileCoefficients>::initializeDirectSolver(
     const bool DirBC_Interior, const StencilDistributionMethod stencil_distribution_method)
 {
-    if (stencil_distribution_method == StencilDistributionMethod::CPU_TAKE) {
+    if (stencil_distribution_method == StencilDistributionMethod::TAKE) {
         op_directSolver_ = std::make_unique<DirectSolverTake<LevelCacheType>>(*grid_, *level_cache_, DirBC_Interior);
     }
-    else if (stencil_distribution_method == StencilDistributionMethod::CPU_GIVE) {
+    else if (stencil_distribution_method == StencilDistributionMethod::GIVE) {
         op_directSolver_ = std::make_unique<DirectSolverGive<LevelCacheType>>(*grid_, *level_cache_, DirBC_Interior);
     }
 
@@ -152,10 +152,10 @@ template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoeff
 void Level<DomainGeometry, DensityProfileCoefficients>::initializeSmoothing(
     const bool DirBC_Interior, const StencilDistributionMethod stencil_distribution_method)
 {
-    if (stencil_distribution_method == StencilDistributionMethod::CPU_TAKE) {
+    if (stencil_distribution_method == StencilDistributionMethod::TAKE) {
         op_smoother_ = std::make_unique<SmootherTake<LevelCacheType>>(*grid_, *level_cache_, DirBC_Interior);
     }
-    else if (stencil_distribution_method == StencilDistributionMethod::CPU_GIVE) {
+    else if (stencil_distribution_method == StencilDistributionMethod::GIVE) {
         op_smoother_ = std::make_unique<SmootherGive<LevelCacheType>>(*grid_, *level_cache_, DirBC_Interior);
     }
     if (!op_smoother_)
@@ -177,11 +177,11 @@ template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoeff
 void Level<DomainGeometry, DensityProfileCoefficients>::initializeExtrapolatedSmoothing(
     const bool DirBC_Interior, const StencilDistributionMethod stencil_distribution_method)
 {
-    if (stencil_distribution_method == StencilDistributionMethod::CPU_TAKE) {
+    if (stencil_distribution_method == StencilDistributionMethod::TAKE) {
         op_extrapolated_smoother_ =
             std::make_unique<ExtrapolatedSmootherTake<LevelCacheType>>(*grid_, *level_cache_, DirBC_Interior);
     }
-    else if (stencil_distribution_method == StencilDistributionMethod::CPU_GIVE) {
+    else if (stencil_distribution_method == StencilDistributionMethod::GIVE) {
         op_extrapolated_smoother_ =
             std::make_unique<ExtrapolatedSmootherGive<LevelCacheType>>(*grid_, *level_cache_, DirBC_Interior);
     }
diff --git a/src/ConfigParser/config_parser.cpp b/src/ConfigParser/config_parser.cpp
index f1cd7f6af..57b14e9b5 100644
--- a/src/ConfigParser/config_parser.cpp
+++ b/src/ConfigParser/config_parser.cpp
@@ -8,7 +8,7 @@ ConfigParser::ConfigParser()
     parser_.add<int>("paraview", '\0', "Generate ParaView output (0/1).", OPTIONAL, 0);
     parser_.add<int>("DirBC_Interior", '\0', "Interior BC type (0=Across-origin, 1=Dirichlet).", OPTIONAL, 0,
                      cmdline::oneof(0, 1));
-    parser_.add<int>("stencilDistributionMethod", '\0', "Stencil distribution (0=CPU_Take,1=CPU_Give)", OPTIONAL, 0,
+    parser_.add<int>("stencilDistributionMethod", '\0', "Stencil distribution (0=TAKE,1=GIVE)", OPTIONAL, 0,
                      cmdline::oneof(0, 1));
     parser_.add<int>("cacheDensityProfileCoefficients", '\0', "Cache density coefficients (0/1).", OPTIONAL, 1,
                      cmdline::oneof(0, 1));
@@ -97,8 +97,8 @@ bool ConfigParser::parse(int argc, char* argv[])
     paraview_             = parser_.get<int>("paraview") != 0;
     DirBC_Interior_       = parser_.get<int>("DirBC_Interior") != 0;
     const int methodValue = parser_.get<int>("stencilDistributionMethod");
-    if (methodValue == static_cast<int>(StencilDistributionMethod::CPU_TAKE) ||
-        methodValue == static_cast<int>(StencilDistributionMethod::CPU_GIVE)) {
+    if (methodValue == static_cast<int>(StencilDistributionMethod::TAKE) ||
+        methodValue == static_cast<int>(StencilDistributionMethod::GIVE)) {
         stencil_distribution_method_ = static_cast<StencilDistributionMethod>(methodValue);
     }
     else {
diff --git a/src/GMGPolar/gmgpolar.cpp b/src/GMGPolar/gmgpolar.cpp
index d237e7d48..2b3d61458 100644
--- a/src/GMGPolar/gmgpolar.cpp
+++ b/src/GMGPolar/gmgpolar.cpp
@@ -11,7 +11,7 @@ IGMGPolar::IGMGPolar(const PolarGrid& grid)
     , paraview_(false)
     // Numerical method setup
     , DirBC_Interior_(true)
-    , stencil_distribution_method_(StencilDistributionMethod::CPU_GIVE)
+    , stencil_distribution_method_(StencilDistributionMethod::GIVE)
     , cache_density_profile_coefficients_(true)
     , cache_domain_geometry_(false)
     // Multigrid settings
diff --git a/src/convergence_order.cpp b/src/convergence_order.cpp
index 99cec9e37..e63d68b12 100644
--- a/src/convergence_order.cpp
+++ b/src/convergence_order.cpp
@@ -16,7 +16,7 @@ int main(int argc, char* argv[])
     const int verbose   = 0;
     const bool paraview = false;
 
-    const StencilDistributionMethod stencilDistributionMethod = StencilDistributionMethod::CPU_TAKE;
+    const StencilDistributionMethod stencilDistributionMethod = StencilDistributionMethod::TAKE;
     const bool cacheDensityProfileCoefficients                = true;
     const bool cacheDomainGeometry                            = true;
 
diff --git a/src/strong_scaling.cpp b/src/strong_scaling.cpp
index fbe390d39..f5ed1f7ef 100644
--- a/src/strong_scaling.cpp
+++ b/src/strong_scaling.cpp
@@ -56,7 +56,7 @@ void runTest(int divideBy2, std::ostream& outfile)
     const int verbose   = 1;
     const bool paraview = false;
 
-    const StencilDistributionMethod stencilDistributionMethod = StencilDistributionMethod::CPU_GIVE;
+    const StencilDistributionMethod stencilDistributionMethod = StencilDistributionMethod::GIVE;
     const bool cacheDensityProfileCoefficients                = true;
     const bool cacheDomainGeometry                            = false;
 
@@ -121,10 +121,10 @@ void runTest(int divideBy2, std::ostream& outfile)
     solver.printTimings();
 
     std::string stencil_string = "";
-    if (solver.stencilDistributionMethod() == StencilDistributionMethod::CPU_TAKE) {
+    if (solver.stencilDistributionMethod() == StencilDistributionMethod::TAKE) {
         stencil_string = "Take";
     }
-    else if (solver.stencilDistributionMethod() == StencilDistributionMethod::CPU_GIVE) {
+    else if (solver.stencilDistributionMethod() == StencilDistributionMethod::GIVE) {
         stencil_string = "Give";
     }
 
diff --git a/src/weak_scaling.cpp b/src/weak_scaling.cpp
index 52279cfeb..9f07e6044 100644
--- a/src/weak_scaling.cpp
+++ b/src/weak_scaling.cpp
@@ -56,7 +56,7 @@ void runTest(int divideBy2, std::ostream& outfile)
     const int verbose   = 1;
     const bool paraview = false;
 
-    const StencilDistributionMethod stencilDistributionMethod = StencilDistributionMethod::CPU_GIVE;
+    const StencilDistributionMethod stencilDistributionMethod = StencilDistributionMethod::GIVE;
     const bool cacheDensityProfileCoefficients                = true;
     const bool cacheDomainGeometry                            = false;
 
@@ -119,10 +119,10 @@ void runTest(int divideBy2, std::ostream& outfile)
     solver.solve(boundary_conditions, source_term);
 
     std::string stencil_string = "";
-    if (solver.stencilDistributionMethod() == StencilDistributionMethod::CPU_TAKE) {
+    if (solver.stencilDistributionMethod() == StencilDistributionMethod::TAKE) {
         stencil_string = "Take";
     }
-    else if (solver.stencilDistributionMethod() == StencilDistributionMethod::CPU_GIVE) {
+    else if (solver.stencilDistributionMethod() == StencilDistributionMethod::GIVE) {
         stencil_string = "Give";
     }
 
diff --git a/tests/GMGPolar/convergence_order.cpp b/tests/GMGPolar/convergence_order.cpp
index b9e8c3ea2..c2c6a77f6 100644
--- a/tests/GMGPolar/convergence_order.cpp
+++ b/tests/GMGPolar/convergence_order.cpp
@@ -173,7 +173,7 @@ get_gmgpolar_error(PolarGrid const& grid, CzarnyGeometry const& domain_geometry,
 
     // --- Discretization and method settings --- //
     gmgpolar.DirBC_Interior(false); // Use across-origin calculation
-    gmgpolar.stencilDistributionMethod(StencilDistributionMethod::CPU_TAKE);
+    gmgpolar.stencilDistributionMethod(StencilDistributionMethod::TAKE);
     gmgpolar.cacheDensityProfileCoefficients(true);
     gmgpolar.cacheDomainGeometry(true);
 
diff --git a/tests/GMGPolar/pcg_tests.cpp b/tests/GMGPolar/pcg_tests.cpp
index 2c3a0ccfa..df7becd50 100644
--- a/tests/GMGPolar/pcg_tests.cpp
+++ b/tests/GMGPolar/pcg_tests.cpp
@@ -156,7 +156,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 1>, // divideBy2
         std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, false>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
@@ -198,7 +198,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, false>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
@@ -240,7 +240,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
@@ -282,7 +282,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_TAKE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::TAKE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
@@ -322,7 +322,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_TAKE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::TAKE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::NONE>, // extrapolation
@@ -362,7 +362,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, false>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
@@ -402,7 +402,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::COMBINED>, // extrapolation
@@ -442,7 +442,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_FULL_GRID_SMOOTHING>, // extrapolation
@@ -481,7 +481,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::NONE>, // extrapolation
@@ -520,7 +520,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_TAKE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::TAKE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
@@ -559,7 +559,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, false>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::NONE>, // extrapolation
@@ -598,7 +598,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 1>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, false>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
diff --git a/tests/GMGPolar/solve_tests.cpp b/tests/GMGPolar/solve_tests.cpp
index 074bb3405..9d1460ba5 100644
--- a/tests/GMGPolar/solve_tests.cpp
+++ b/tests/GMGPolar/solve_tests.cpp
@@ -138,7 +138,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 1>, // divideBy2
         std::integral_constant<int, 1>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, false>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
@@ -174,7 +174,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, false>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
@@ -210,7 +210,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
@@ -246,7 +246,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_TAKE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::TAKE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
@@ -280,7 +280,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_TAKE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::TAKE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::NONE>, // extrapolation
@@ -314,7 +314,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, false>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
@@ -348,7 +348,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::COMBINED>, // extrapolation
@@ -382,7 +382,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, false>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_FULL_GRID_SMOOTHING>, // extrapolation
@@ -415,7 +415,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, true>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::NONE>, // extrapolation
@@ -448,7 +448,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_TAKE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::TAKE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, true>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
@@ -481,7 +481,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, false>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::NONE>, // extrapolation
@@ -514,7 +514,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 1>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, false>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::IMPLICIT_EXTRAPOLATION>, // extrapolation
@@ -547,7 +547,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, false>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::NONE>, // extrapolation
@@ -580,7 +580,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, false>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::NONE>, // extrapolation
@@ -613,7 +613,7 @@ using gmgpolar_test_suite = testing::Types<
         std::integral_constant<int, 0>, // divideBy2
         std::integral_constant<int, 0>, // verbose
         std::integral_constant<bool, false>, // DirBC_Interior
-        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::CPU_GIVE>, // StencilDistributionMethod
+        std::integral_constant<StencilDistributionMethod, StencilDistributionMethod::GIVE>, // StencilDistributionMethod
         std::integral_constant<bool, true>, // cacheDensityProfileCoefficient
         std::integral_constant<bool, false>, // cacheDomainGeometry
         std::integral_constant<ExtrapolationType, ExtrapolationType::NONE>, // extrapolation