WIP: More porting

include/DirectSolver/DirectSolverGive/directSolverGive.h

@@ -13,7 +13,7 @@ class DirectSolverGive : public DirectSolver<LevelCacheType>
                               bool DirBC_Interior);
 
     // Note: The rhs (right-hand side) vector gets overwritten during the solution process.
-    void solveInPlace(HostVector<double> solution) override;
+    void solveInPlace(Vector<double> solution) override;
 
 private:
 #ifdef GMGPOLAR_USE_MUMPS

include/DirectSolver/DirectSolverGive/directSolverGive.inl

@@ -14,9 +14,8 @@ DirectSolverGive<LevelCacheType>::DirectSolverGive(const PolarGrid<DefaultMemory
 }
 
 template <class LevelCacheType>
-void DirectSolverGive<LevelCacheType>::solveInPlace(HostVector<double> h_solution)
+void DirectSolverGive<LevelCacheType>::solveInPlace(Vector<double> solution)
 {
-    auto solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_solution);
     // Adjusts the right-hand side vector to account for symmetry corrections.
     // This transforms the system matrixA * solution = rhs into the equivalent system:
     // symmetric_DBc(matrixA) * solution = rhs - applySymmetryShift(rhs).
@@ -25,6 +24,4 @@ void DirectSolverGive<LevelCacheType>::solveInPlace(HostVector<double> h_solutio
     applySymmetryShift(solution);
     // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver.
     system_solver_.solveInPlace(solution);
-
-    Kokkos::deep_copy(h_solution, solution);
 }

include/DirectSolver/DirectSolverTake/directSolverTake.h

@@ -13,7 +13,7 @@ class DirectSolverTake : public DirectSolver<LevelCacheType>
                               bool DirBC_Interior);
 
     // Note: The rhs (right-hand side) vector gets overwritten during the solution process.
-    void solveInPlace(HostVector<double> solution) override;
+    void solveInPlace(Vector<double> solution) override;
 
 private:
 #ifdef GMGPOLAR_USE_MUMPS

include/DirectSolver/DirectSolverTake/directSolverTake.inl

@@ -14,9 +14,8 @@ DirectSolverTake<LevelCacheType>::DirectSolverTake(const PolarGrid<DefaultMemory
 }
 
 template <class LevelCacheType>
-void DirectSolverTake<LevelCacheType>::solveInPlace(HostVector<double> h_solution)
+void DirectSolverTake<LevelCacheType>::solveInPlace(Vector<double> solution)
 {
-    auto solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_solution);
     // Adjusts the right-hand side vector to account for symmetry corrections.
     // This transforms the system matrixA * solution = rhs into the equivalent system:
     // symmetric_DBc(matrixA) * solution = rhs - applySymmetryShift(rhs).
@@ -25,6 +24,4 @@ void DirectSolverTake<LevelCacheType>::solveInPlace(HostVector<double> h_solutio
     applySymmetryShift(solution);
     // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver.
     system_solver_.solveInPlace(solution);
-
-    Kokkos::deep_copy(h_solution, solution);
 }

include/DirectSolver/directSolver.h

@@ -37,7 +37,7 @@ class DirectSolver
     virtual ~DirectSolver() = default;
 
     // Note: The rhs (right-hand side) vector gets overwritten during the solution process.
-    virtual void solveInPlace(HostVector<double> solution) = 0;
+    virtual void solveInPlace(Vector<double> solution) = 0;
 
 protected:
     const PolarGrid<DefaultMemorySpace>& grid_;

include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h

@@ -68,7 +68,7 @@ class ExtrapolatedSmootherGive : public ExtrapolatedSmoother<LevelCacheType>
     // Parallel implementation using OpenMP:
     // Scedule every 2nd/4th line in parallel to avoid race conditions arising from the A-Give distribution.
     // Sceduling every 3rd line in parallel would also be possible, but is less natural for the 2 coloring.
-    void extrapolatedSmoothing(HostVector<double> x, HostConstVector<double> rhs, HostVector<double> temp) override;
+    void extrapolatedSmoothing(Vector<double> x, ConstVector<double> rhs, Vector<double> temp) override;
 
 private:
     /* ------------------- */

include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.inl

@@ -42,14 +42,9 @@ ExtrapolatedSmootherGive<LevelCacheType>::ExtrapolatedSmootherGive(const PolarGr
 //   - The system is then solved in-place in temp, and the results
 //     are copied back to x.
 template <class LevelCacheType>
-void ExtrapolatedSmootherGive<LevelCacheType>::extrapolatedSmoothing(HostVector<double> h_x,
-                                                                     HostConstVector<double> h_rhs,
-                                                                     HostVector<double> h_temp)
+void ExtrapolatedSmootherGive<LevelCacheType>::extrapolatedSmoothing(Vector<double> x, ConstVector<double> rhs,
+                                                                     Vector<double> temp)
 {
-    auto x    = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_x);
-    auto rhs  = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs);
-    auto temp = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_temp);
-
     assert(x.size() == rhs.size());
     assert(temp.size() == rhs.size());
 
@@ -91,7 +86,4 @@ void ExtrapolatedSmootherGive<LevelCacheType>::extrapolatedSmoothing(HostVector<
 
     applyAscOrthoWhiteRadialSection(x, rhs, temp);
     solveWhiteRadialSection(x, temp);
-
-    Kokkos::deep_copy(h_x, x);
-    Kokkos::deep_copy(h_temp, temp);
 }

include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.h

@@ -66,7 +66,7 @@ class ExtrapolatedSmootherTake : public ExtrapolatedSmoother<LevelCacheType>
     // Performs one full coupled extrapolated smoothing sweep:
     //   BC -> WC -> BR -> WR
     // using temp as RHS workspace.
-    void extrapolatedSmoothing(HostVector<double> x, HostConstVector<double> rhs, HostVector<double> temp) override;
+    void extrapolatedSmoothing(Vector<double> x, ConstVector<double> rhs, Vector<double> temp) override;
 
 private:
     /* ------------------- */

include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/extrapolatedSmootherTake.inl

@@ -43,14 +43,9 @@ ExtrapolatedSmootherTake<LevelCacheType>::ExtrapolatedSmootherTake(const PolarGr
 //     are copied back to x.
 
 template <class LevelCacheType>
-void ExtrapolatedSmootherTake<LevelCacheType>::extrapolatedSmoothing(HostVector<double> h_x,
-                                                                     HostConstVector<double> h_rhs,
-                                                                     HostVector<double> h_temp)
+void ExtrapolatedSmootherTake<LevelCacheType>::extrapolatedSmoothing(Vector<double> x, ConstVector<double> rhs,
+                                                                     Vector<double> temp)
 {
-    auto x    = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_x);
-    auto rhs  = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs);
-    auto temp = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_temp);
-
     assert(x.size() == rhs.size());
     assert(temp.size() == rhs.size());
 
@@ -65,7 +60,4 @@ void ExtrapolatedSmootherTake<LevelCacheType>::extrapolatedSmoothing(HostVector<
 
     applyAscOrthoWhiteRadialSection(x, rhs, temp);
     solveWhiteRadialSection(x, temp);
-
-    Kokkos::deep_copy(h_x, x);
-    Kokkos::deep_copy(h_temp, temp);
 }

include/ExtrapolatedSmoother/extrapolatedSmoother.h

@@ -36,7 +36,7 @@ class ExtrapolatedSmoother
     }
     virtual ~ExtrapolatedSmoother() = default;
 
-    virtual void extrapolatedSmoothing(HostVector<double> x, HostConstVector<double> rhs, HostVector<double> temp) = 0;
+    virtual void extrapolatedSmoothing(Vector<double> x, ConstVector<double> rhs, Vector<double> temp) = 0;
 
 protected:
     const PolarGrid<DefaultMemorySpace>& grid_;

include/GMGPolar/MultigridMethods/extrapolated_multigrid_F_Cycle.h

@@ -2,9 +2,9 @@
 
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
 void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigrid_F_Cycle(int level_depth,
-                                                                                          HostVector<double> solution,
-                                                                                          HostConstVector<double> rhs,
-                                                                                          HostVector<double> residual)
+                                                                                          HostVector<double> h_solution,
+                                                                                          HostConstVector<double> h_rhs,
+                                                                                          HostVector<double> h_residual)
 {
     assert(level_depth == 0);
     assert(number_of_levels_ >= 2);
@@ -25,6 +25,12 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
     /* ------------ */
     auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
 
+    auto residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_residual);
+    auto solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_solution);
+    auto rhs      = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs); // const
+    auto next_level_residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.residual());
+    auto next_level_solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.solution());
+
     for (int i = 0; i < pre_smoothing_steps_; i++) {
         if (level.level_depth() == 0 && !full_grid_smoothing_) {
             level.extrapolatedSmoothing(solution, rhs, residual);
@@ -38,21 +44,23 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
     t_avg_MGC_preSmoothing_ += std::chrono::duration<double>(end_MGC_preSmoothing - start_MGC_preSmoothing).count();
 
     /* -------------------------------------------- */
-    /* Compute the restricted extrapolated residual */
+    /* Compute the restricted extrapolated h_residual */
     /* -------------------------------------------- */
     auto start_MGC_residual = std::chrono::high_resolution_clock::now();
 
     // P_ex^T (f_l - A_l*u_l)
     level.computeResidual(residual, rhs, solution);
-    extrapolatedRestriction(level.level_depth(), next_level.residual(), residual);
+    extrapolatedRestriction(level.level_depth(), next_level_residual, residual);
 
     // f_{l-1} - A_{l-1}* Inject(u_l)
-    injection(level.level_depth(), next_level.solution(), solution);
-    next_level.computeResidual(next_level.error_correction(), next_level.rhs(), next_level.solution());
+    injection(level.level_depth(), next_level_solution, solution);
+    next_level.computeResidual(next_level.error_correction(), next_level.rhs(), next_level_solution);
 
     // res_ex = 4/3 * P_ex^T (f_l - A_l*u_l) - 1/3 * (f_{l-1} - A_{l-1}* Inject(u_l))
-    linear_combination(next_level.residual(), 4.0 / 3.0, HostConstVector<double>(next_level.error_correction()),
+    linear_combination(next_level_residual, 4.0 / 3.0, ConstVector<double>(next_level.error_correction()),
                        -1.0 / 3.0);
+    Kokkos::deep_copy(next_level.residual(), next_level_residual);
+    Kokkos::deep_copy(next_level.solution(), next_level_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();
@@ -66,12 +74,13 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
     assign(next_level.error_correction(), 0.0);
 
     /* Solve for the error by recursively calling the multigrid cycle. */
-    multigrid_F_Cycle(next_level.level_depth(), next_level.error_correction(), next_level.residual(),
+	auto h_next_level_error_correction = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, next_level.error_correction());
+    multigrid_F_Cycle(next_level.level_depth(), h_next_level_error_correction, next_level.residual(),
                       next_level.solution());
 
     /* Don't do a second recursion on the coarsest level since the DirectSolver is exact. */
     if (next_level.level_depth() != number_of_levels_ - 1) {
-        multigrid_V_Cycle(next_level.level_depth(), next_level.error_correction(), next_level.residual(),
+        multigrid_V_Cycle(next_level.level_depth(), h_next_level_error_correction, next_level.residual(),
                           next_level.solution());
     }
 
@@ -80,12 +89,13 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
     /* -------------------------- */
     // Use 'residual' instead of 'level.error_correction()' as a temporary buffer.
     // Note: 'level.error_correction()' has size 0 at level depth = 0.
-    extrapolatedProlongation(next_level.level_depth(), residual, next_level.error_correction());
+    Kokkos::deep_copy(next_level.error_correction(), h_next_level_error_correction);
+	extrapolatedProlongation(next_level.level_depth(), residual, next_level.error_correction());
 
     /* ----------------------------------- */
     /* Compute the corrected approximation */
     /* ----------------------------------- */
-    add(solution, HostConstVector<double>(residual));
+    add(solution, ConstVector<double>(residual));
 
     /* ------------- */
     /* Postsmoothing */
@@ -100,6 +110,8 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
             level.smoothing(solution, rhs, residual);
         }
     }
+	Kokkos::deep_copy(h_residual, residual);
+	Kokkos::deep_copy(h_solution, solution);
 
     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();

include/GMGPolar/MultigridMethods/extrapolated_multigrid_V_Cycle.h

@@ -2,9 +2,9 @@
 
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
 void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigrid_V_Cycle(int level_depth,
-                                                                                          HostVector<double> solution,
-                                                                                          HostConstVector<double> rhs,
-                                                                                          HostVector<double> residual)
+                                                                                          HostVector<double> h_solution,
+                                                                                          HostConstVector<double> h_rhs,
+                                                                                          HostVector<double> h_residual)
 {
     assert(level_depth == 0);
     assert(number_of_levels_ >= 2);
@@ -25,6 +25,12 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
     /* ------------ */
     auto start_MGC_preSmoothing = std::chrono::high_resolution_clock::now();
 
+    auto residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_residual);
+    auto solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_solution);
+    auto rhs      = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), h_rhs); // const
+    auto next_level_residual = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.residual());
+    auto next_level_solution = Kokkos::create_mirror_view_and_copy(DefaultMemorySpace(), next_level.solution());
+
     for (int i = 0; i < pre_smoothing_steps_; i++) {
         if (level.level_depth() == 0 && !full_grid_smoothing_) {
             level.extrapolatedSmoothing(solution, rhs, residual);
@@ -44,15 +50,17 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
 
     // P_ex^T (f_l - A_l*u_l)
     level.computeResidual(residual, rhs, solution);
-    extrapolatedRestriction(level.level_depth(), next_level.residual(), residual);
+    extrapolatedRestriction(level.level_depth(), next_level_residual, residual);
 
     // f_{l-1} - A_{l-1}* Inject(u_l)
-    injection(level.level_depth(), next_level.solution(), solution);
-    next_level.computeResidual(next_level.error_correction(), next_level.rhs(), next_level.solution());
+    injection(level.level_depth(), next_level_solution, solution);
+    next_level.computeResidual(next_level.error_correction(), next_level.rhs(), next_level_solution);
 
     // res_ex = 4/3 * P_ex^T (f_l - A_l*u_l) - 1/3 * (f_{l-1} - A_{l-1}* Inject(u_l))
-    linear_combination(next_level.residual(), 4.0 / 3.0, HostConstVector<double>(next_level.error_correction()),
+    linear_combination(next_level_residual, 4.0 / 3.0, ConstVector<double>(next_level.error_correction()),
                        -1.0 / 3.0);
+    Kokkos::deep_copy(next_level.residual(), next_level_residual);
+    Kokkos::deep_copy(next_level.solution(), next_level_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();
@@ -66,20 +74,22 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
     assign(next_level.error_correction(), 0.0);
 
     /* Solve for the error by recursively calling the multigrid cycle. */
-    multigrid_V_Cycle(next_level.level_depth(), next_level.error_correction(), next_level.residual(),
+	auto h_next_level_error_correction = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, next_level.error_correction());
+    multigrid_V_Cycle(next_level.level_depth(), h_next_level_error_correction, next_level.residual(),
                       next_level.solution());
 
     /* -------------------------- */
     /* Interpolate the correction */
     /* -------------------------- */
     // Use 'residual' instead of 'level.error_correction()' as a temporary buffer.
     // Note: 'level.error_correction()' has size 0 at level depth = 0.
+	Kokkos::deep_copy(next_level.error_correction(), h_next_level_error_correction);
     extrapolatedProlongation(next_level.level_depth(), residual, next_level.error_correction());
 
     /* ----------------------------------- */
     /* Compute the corrected approximation */
     /* ----------------------------------- */
-    add(solution, HostConstVector<double>(residual));
+    add(solution, ConstVector<double>(residual));
 
     /* ------------- */
     /* Postsmoothing */
@@ -94,6 +104,8 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::extrapolated_multigri
             level.smoothing(solution, rhs, residual);
         }
     }
+	Kokkos::deep_copy(h_residual, residual);
+	Kokkos::deep_copy(h_solution, solution);
 
     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();