WIP: Start GPU porting

include/ConfigParser/config_parser.h

@@ -46,7 +46,7 @@ class ConfigParser
     bool cacheDensityProfileCoefficients() const;
     bool cacheDomainGeometry() const;
 
-    const PolarGrid& grid() const;
+    const PolarGrid<Kokkos::HostSpace>& grid() const;
 
     // Full Multigrid Method
     bool FMG() const;
@@ -96,7 +96,7 @@ class ConfigParser
     bool cache_density_profile_coefficients_;
     bool cache_domain_geometry_;
     // Grid configuration
-    PolarGrid grid_;
+    PolarGrid<Kokkos::HostSpace> grid_;
     // Multigrid settings
     ExtrapolationType extrapolation_;
     int max_levels_;

include/DirectSolver/DirectSolverGive/applySymmetryShift.inl

@@ -5,16 +5,16 @@
 /* ----------------------- */
 
 template <class LevelCacheType>
-void DirectSolverGive<LevelCacheType>::applySymmetryShiftInnerBoundary(HostVector<double> x) const
+void DirectSolverGive<LevelCacheType>::applySymmetryShiftInnerBoundary(Vector<double> x) const
 {
-    const PolarGrid& grid             = DirectSolver<LevelCacheType>::grid_;
-    const LevelCacheType& level_cache = DirectSolver<LevelCacheType>::level_cache_;
+    const PolarGrid<DefaultMemorySpace>& grid = DirectSolver<LevelCacheType>::grid_;
+    const LevelCacheType& level_cache         = DirectSolver<LevelCacheType>::level_cache_;
 
     assert(DirectSolver<LevelCacheType>::DirBC_Interior_);
 
     Kokkos::parallel_for(
-        "DirectSolverGive: applySymmetryShiftInnerBoundary",
-        Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, 1), KOKKOS_LAMBDA(const int) {
+        "DirectSolverGive: applySymmetryShiftInnerBoundary", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, 1),
+        KOKKOS_LAMBDA(const int) {
             int i_r;
             double r;
             int global_index;
@@ -70,14 +70,14 @@ void DirectSolverGive<LevelCacheType>::applySymmetryShiftInnerBoundary(HostVecto
 }
 
 template <class LevelCacheType>
-void DirectSolverGive<LevelCacheType>::applySymmetryShiftOuterBoundary(HostVector<double> x) const
+void DirectSolverGive<LevelCacheType>::applySymmetryShiftOuterBoundary(Vector<double> x) const
 {
-    const PolarGrid& grid             = DirectSolver<LevelCacheType>::grid_;
-    const LevelCacheType& level_cache = DirectSolver<LevelCacheType>::level_cache_;
+    const PolarGrid<DefaultMemorySpace>& grid = DirectSolver<LevelCacheType>::grid_;
+    const LevelCacheType& level_cache         = DirectSolver<LevelCacheType>::level_cache_;
 
     Kokkos::parallel_for(
-        "DirectSolverGive: applySymmetryShiftOuterBoundary",
-        Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, 1), KOKKOS_LAMBDA(const int) {
+        "DirectSolverGive: applySymmetryShiftOuterBoundary", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, 1),
+        KOKKOS_LAMBDA(const int) {
             int i_r;
             double r;
             int global_index;
@@ -133,10 +133,10 @@ void DirectSolverGive<LevelCacheType>::applySymmetryShiftOuterBoundary(HostVecto
 }
 
 template <class LevelCacheType>
-void DirectSolverGive<LevelCacheType>::applySymmetryShift(HostVector<double> x) const
+void DirectSolverGive<LevelCacheType>::applySymmetryShift(Vector<double> x) const
 {
-    const PolarGrid& grid     = DirectSolver<LevelCacheType>::grid_;
-    const bool DirBC_Interior = DirectSolver<LevelCacheType>::DirBC_Interior_;
+    const PolarGrid<DefaultMemorySpace>& grid = DirectSolver<LevelCacheType>::grid_;
+    const bool DirBC_Interior                 = DirectSolver<LevelCacheType>::DirBC_Interior_;
 
     assert(std::ssize(x) == grid.numberOfNodes());
     assert(grid.nr() >= 4);

include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl

@@ -7,19 +7,17 @@ namespace direct_solver_give
 
 #ifdef GMGPOLAR_USE_MUMPS
 // When using the MUMPS solver, the matrix is assembled in COO format.
-static KOKKOS_INLINE_FUNCTION void updateMatrixElement(const SparseMatrixCOO<double, Kokkos::HostSpace>& matrix,
-                                                       int ptr, const int offset, const int row, const int column,
-                                                       const double value)
+static KOKKOS_INLINE_FUNCTION void updateMatrixElement(const SparseMatrixCOO<double>& matrix, int ptr, const int offset,
+                                                       const int row, const int column, const double value)
 {
     matrix.set_row_index(ptr + offset, row);
     matrix.set_col_index(ptr + offset, column);
     matrix.increase_value(ptr + offset, value);
 }
 #else
 // When using the in-house solver, the matrix is stored in CSR format.
-static KOKKOS_INLINE_FUNCTION void updateMatrixElement(const SparseMatrixCSR<double, Kokkos::HostSpace>& matrix,
-                                                       int ptr, const int offset, const int row, const int column,
-                                                       const double value)
+static KOKKOS_INLINE_FUNCTION void updateMatrixElement(const SparseMatrixCSR<double>& matrix, int ptr, const int offset,
+                                                       const int row, const int column, const double value)
 {
     matrix.set_row_nz_index(row, offset, column);
     matrix.increase_row_nz_entry(row, offset, value);
@@ -28,8 +26,9 @@ static KOKKOS_INLINE_FUNCTION void updateMatrixElement(const SparseMatrixCSR<dou
 
 template <typename LevelCacheType, typename SystemMatrix>
 static KOKKOS_INLINE_FUNCTION void
-nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& grid, const LevelCacheType& level_cache,
-                          const bool DirBC_Interior, const SystemMatrix& solver_matrix)
+nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid<DefaultMemorySpace>& grid,
+                          const LevelCacheType& level_cache, const bool DirBC_Interior,
+                          const SystemMatrix& solver_matrix)
 {
     /* ---------------------------------------- */
     /* Compute or retrieve stencil coefficients */
@@ -799,24 +798,24 @@ typename DirectSolverGive<LevelCacheType>::SystemMatrix DirectSolverGive<LevelCa
     using direct_solver_give::nodeBuildSolverMatrixGive;
     using direct_solver_give::validateSolverMatrixIndexing;
 
-    const PolarGrid& grid             = DirectSolver<LevelCacheType>::grid_;
-    const LevelCacheType& level_cache = DirectSolver<LevelCacheType>::level_cache_;
-    const bool DirBC_Interior         = DirectSolver<LevelCacheType>::DirBC_Interior_;
+    const PolarGrid<DefaultMemorySpace>& grid = DirectSolver<LevelCacheType>::grid_;
+    const LevelCacheType& level_cache         = DirectSolver<LevelCacheType>::level_cache_;
+    const bool DirBC_Interior                 = DirectSolver<LevelCacheType>::DirBC_Interior_;
 
     assert(validateSolverMatrixIndexing(grid, DirBC_Interior) && "Solver matrix indexing is inconsistent");
 
     const int n = grid.numberOfNodes();
 
 #ifdef GMGPOLAR_USE_MUMPS
     const int nnz = getNonZeroCountSolverMatrix(grid, DirBC_Interior);
-    SparseMatrixCOO<double, Kokkos::HostSpace> solver_matrix(n, n, nnz);
+    SparseMatrixCOO<double> solver_matrix(n, n, nnz);
     solver_matrix.is_symmetric(true);
 #else
     std::function<int(int)> nnz_per_row = [&](int global_index) {
         return getStencilSize(global_index, grid, DirBC_Interior);
     };
 
-    SparseMatrixCSR<double, Kokkos::HostSpace> solver_matrix(n, n, nnz_per_row);
+    SparseMatrixCSR<double> solver_matrix(n, n, nnz_per_row);
 #endif
 
     /* ---------------- */
@@ -830,7 +829,7 @@ typename DirectSolverGive<LevelCacheType>::SystemMatrix DirectSolverGive<LevelCa
         const int num_circular_tasks = (num_circle_tasks - start_circle + 2) / 3;
         Kokkos::parallel_for(
             "DirectSolverGive: BuildSolverMatrix (Circular)",
-            Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, num_circular_tasks),
+            Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, num_circular_tasks),
             KOKKOS_LAMBDA(const int circle_task) {
                 const int i_r = start_circle + circle_task * 3;
                 for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) {
@@ -853,7 +852,7 @@ typename DirectSolverGive<LevelCacheType>::SystemMatrix DirectSolverGive<LevelCa
     for (int i_theta = 0; i_theta < additional_radial_tasks; i_theta++) {
         Kokkos::parallel_for(
             "DirectSolverGive: BuildSolverMatrix (Radial, additional)",
-            Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, 1), KOKKOS_LAMBDA(const int) {
+            Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, 1), KOKKOS_LAMBDA(const int) {
                 for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) {
                     nodeBuildSolverMatrixGive(i_r, i_theta, grid, level_cache, DirBC_Interior, solver_matrix);
                 }
@@ -865,7 +864,7 @@ typename DirectSolverGive<LevelCacheType>::SystemMatrix DirectSolverGive<LevelCa
         const int num_radial_batches = (num_radial_tasks - start_radial + 2) / 3;
         Kokkos::parallel_for(
             "DirectSolverGive: BuildSolverMatrix (Radial)",
-            Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, num_radial_batches),
+            Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, num_radial_batches),
             KOKKOS_LAMBDA(const int radial_task) {
                 const int i_theta = additional_radial_tasks + start_radial + radial_task * 3;
                 for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) {

include/DirectSolver/DirectSolverGive/directSolverGive.h

@@ -9,18 +9,19 @@ template <class LevelCacheType>
 class DirectSolverGive : public DirectSolver<LevelCacheType>
 {
 public:
-    explicit DirectSolverGive(const PolarGrid& grid, const LevelCacheType& level_cache, bool DirBC_Interior);
+    explicit DirectSolverGive(const PolarGrid<DefaultMemorySpace>& grid, const LevelCacheType& level_cache,
+                              bool DirBC_Interior);
 
     // Note: The rhs (right-hand side) vector gets overwritten during the solution process.
     void solveInPlace(HostVector<double> solution) override;
 
 private:
 #ifdef GMGPOLAR_USE_MUMPS
-    using SystemMatrix = SparseMatrixCOO<double, Kokkos::HostSpace>;
+    using SystemMatrix = SparseMatrixCOO<double>;
     using SystemSolver = CooMumpsSolver;
 #else
-    using SystemMatrix = SparseMatrixCSR<double, Kokkos::HostSpace>;
-    using SystemSolver = SparseLUSolver<double, Kokkos::HostSpace>;
+    using SystemMatrix = SparseMatrixCSR<double>;
+    using SystemSolver = SparseLUSolver<double>;
     // Stored only for the in-house solver (CSR).
     SystemMatrix system_matrix_;
 #endif
@@ -39,9 +40,9 @@ class DirectSolverGive : public DirectSolver<LevelCacheType>
     //    symmetric_DBc(A) * solution = rhs - applySymmetryShift(rhs).
     // The correction modifies the rhs to account for the influence of the Dirichlet boundary conditions,
     // ensuring that the solution at the boundary is correctly adjusted and maintains the required symmetry.
-    void applySymmetryShift(HostVector<double> rhs) const;
-    void applySymmetryShiftInnerBoundary(HostVector<double> x) const;
-    void applySymmetryShiftOuterBoundary(HostVector<double> x) const;
+    void applySymmetryShift(Vector<double> rhs) const;
+    void applySymmetryShiftInnerBoundary(Vector<double> x) const;
+    void applySymmetryShiftOuterBoundary(Vector<double> x) const;
 };
 
 #include "applySymmetryShift.inl"

include/DirectSolver/DirectSolverGive/directSolverGive.inl

@@ -1,8 +1,8 @@
 #pragma once
 
 template <class LevelCacheType>
-DirectSolverGive<LevelCacheType>::DirectSolverGive(const PolarGrid& grid, const LevelCacheType& level_cache,
-                                                   bool DirBC_Interior)
+DirectSolverGive<LevelCacheType>::DirectSolverGive(const PolarGrid<DefaultMemorySpace>& grid,
+                                                   const LevelCacheType& level_cache, bool DirBC_Interior)
     : DirectSolver<LevelCacheType>(grid, level_cache, DirBC_Interior)
 #ifdef GMGPOLAR_USE_MUMPS
     , system_solver_(buildSolverMatrix())
@@ -14,8 +14,9 @@ DirectSolverGive<LevelCacheType>::DirectSolverGive(const PolarGrid& grid, const
 }
 
 template <class LevelCacheType>
-void DirectSolverGive<LevelCacheType>::solveInPlace(HostVector<double> solution)
+void DirectSolverGive<LevelCacheType>::solveInPlace(HostVector<double> h_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).
@@ -24,4 +25,6 @@ void DirectSolverGive<LevelCacheType>::solveInPlace(HostVector<double> solution)
     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/DirectSolverGive/matrixStencil.inl

@@ -3,7 +3,8 @@
 namespace direct_solver_give
 {
 
-static KOKKOS_INLINE_FUNCTION int getStencilSize(int global_index, const PolarGrid& grid, const bool DirBC_Interior)
+static KOKKOS_INLINE_FUNCTION int getStencilSize(int global_index, const PolarGrid<DefaultMemorySpace>& grid,
+                                                 const bool DirBC_Interior)
 {
     int i_r, i_theta;
     grid.multiIndex(global_index, i_r, i_theta);
@@ -33,7 +34,8 @@ static KOKKOS_INLINE_FUNCTION int getStencilSize(int global_index, const PolarGr
     Kokkos::abort("Invalid index for stencil");
 }
 
-static KOKKOS_INLINE_FUNCTION const Stencil& getStencil(const int i_r, const PolarGrid& grid, const bool DirBC_Interior)
+static KOKKOS_INLINE_FUNCTION const Stencil& getStencil(const int i_r, const PolarGrid<DefaultMemorySpace>& grid,
+                                                        const bool DirBC_Interior)
 {
     assert(0 <= i_r && i_r < grid.nr());
     assert(grid.nr() >= 4);
@@ -85,7 +87,8 @@ static KOKKOS_INLINE_FUNCTION const Stencil& getStencil(const int i_r, const Pol
     Kokkos::abort("Invalid index for stencil");
 }
 
-static KOKKOS_INLINE_FUNCTION int getNonZeroCountSolverMatrix(const PolarGrid& grid, const bool DirBC_Interior)
+static KOKKOS_INLINE_FUNCTION int getNonZeroCountSolverMatrix(const PolarGrid<DefaultMemorySpace>& grid,
+                                                              const bool DirBC_Interior)
 {
     const int size_stencil_inner_boundary      = DirBC_Interior ? 1 : 7;
     const int size_stencil_next_inner_boundary = DirBC_Interior ? 6 : 9;
@@ -102,7 +105,8 @@ static KOKKOS_INLINE_FUNCTION int getNonZeroCountSolverMatrix(const PolarGrid& g
 
 /* ----------------------------------------------------------------- */
 /* If the indexing is not smoother-based, please adjust the indexing */
-static KOKKOS_INLINE_FUNCTION int getSolverMatrixIndex(const int i_r, const int i_theta, const PolarGrid& grid,
+static KOKKOS_INLINE_FUNCTION int getSolverMatrixIndex(const int i_r, const int i_theta,
+                                                       const PolarGrid<DefaultMemorySpace>& grid,
                                                        const bool DirBC_Interior)
 {
     const int size_stencil_inner_boundary      = DirBC_Interior ? 1 : 7;
@@ -182,7 +186,8 @@ static KOKKOS_INLINE_FUNCTION int getSolverMatrixIndex(const int i_r, const int
     Kokkos::abort("Invalid index for stencil");
 }
 
-static KOKKOS_INLINE_FUNCTION bool validateSolverMatrixIndexing(const PolarGrid& grid, const bool DirBC_Interior)
+static KOKKOS_INLINE_FUNCTION bool validateSolverMatrixIndexing(const PolarGrid<DefaultMemorySpace>& grid,
+                                                                const bool DirBC_Interior)
 {
     // 1. Check each node: getSolverMatrixIndex == cumulative sum of prior stencil sizes
     for (int global_index = 0; global_index < grid.numberOfNodes(); ++global_index) {

include/DirectSolver/DirectSolverTake/applySymmetryShift.inl

@@ -5,26 +5,26 @@
 /* ----------------------- */
 
 template <class LevelCacheType>
-void DirectSolverTake<LevelCacheType>::applySymmetryShiftInnerBoundary(HostVector<double> x) const
+void DirectSolverTake<LevelCacheType>::applySymmetryShiftInnerBoundary(Vector<double> x) const
 {
-    const PolarGrid& grid             = DirectSolver<LevelCacheType>::grid_;
-    const LevelCacheType& level_cache = DirectSolver<LevelCacheType>::level_cache_;
+    const PolarGrid<DefaultMemorySpace>& grid = DirectSolver<LevelCacheType>::grid_;
+    const LevelCacheType& level_cache         = DirectSolver<LevelCacheType>::level_cache_;
 
     assert(DirectSolver<LevelCacheType>::DirBC_Interior_);
 
     assert(level_cache.cacheDensityProfileCoefficients());
     assert(level_cache.cacheDomainGeometry());
 
-    HostConstVector<double> arr = level_cache.arr();
-    HostConstVector<double> att = level_cache.att();
-    HostConstVector<double> art = level_cache.art();
+    ConstVector<double> arr = level_cache.arr();
+    ConstVector<double> att = level_cache.att();
+    ConstVector<double> art = level_cache.art();
 
     const int i_r    = 1;
     const int ntheta = grid.ntheta();
 
     Kokkos::parallel_for(
         "DirectSolverTake: applySymmetryShiftInnerBoundary",
-        Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, ntheta), KOKKOS_LAMBDA(const int i_theta) {
+        Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, ntheta), KOKKOS_LAMBDA(const int i_theta) {
             const double h1 = grid.radialSpacing(i_r - 1);
             const double k1 = grid.angularSpacing(i_theta - 1);
             const double k2 = grid.angularSpacing(i_theta);
@@ -50,24 +50,24 @@ void DirectSolverTake<LevelCacheType>::applySymmetryShiftInnerBoundary(HostVecto
 }
 
 template <class LevelCacheType>
-void DirectSolverTake<LevelCacheType>::applySymmetryShiftOuterBoundary(HostVector<double> x) const
+void DirectSolverTake<LevelCacheType>::applySymmetryShiftOuterBoundary(Vector<double> x) const
 {
-    const PolarGrid& grid             = DirectSolver<LevelCacheType>::grid_;
-    const LevelCacheType& level_cache = DirectSolver<LevelCacheType>::level_cache_;
+    const PolarGrid<DefaultMemorySpace>& grid = DirectSolver<LevelCacheType>::grid_;
+    const LevelCacheType& level_cache         = DirectSolver<LevelCacheType>::level_cache_;
 
     assert(level_cache.cacheDensityProfileCoefficients());
     assert(level_cache.cacheDomainGeometry());
 
-    HostConstVector<double> arr = level_cache.arr();
-    HostConstVector<double> att = level_cache.att();
-    HostConstVector<double> art = level_cache.art();
+    ConstVector<double> arr = level_cache.arr();
+    ConstVector<double> att = level_cache.att();
+    ConstVector<double> art = level_cache.art();
 
     const int i_r    = grid.nr() - 2;
     const int ntheta = grid.ntheta();
 
     Kokkos::parallel_for(
         "DirectSolverTake: applySymmetryShiftOuterBoundary",
-        Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, ntheta), KOKKOS_LAMBDA(const int i_theta) {
+        Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, ntheta), KOKKOS_LAMBDA(const int i_theta) {
             const double h2 = grid.radialSpacing(i_r);
             const double k1 = grid.angularSpacing(i_theta - 1);
             const double k2 = grid.angularSpacing(i_theta);
@@ -93,10 +93,10 @@ void DirectSolverTake<LevelCacheType>::applySymmetryShiftOuterBoundary(HostVecto
 }
 
 template <class LevelCacheType>
-void DirectSolverTake<LevelCacheType>::applySymmetryShift(HostVector<double> x) const
+void DirectSolverTake<LevelCacheType>::applySymmetryShift(Vector<double> x) const
 {
-    const PolarGrid& grid     = DirectSolver<LevelCacheType>::grid_;
-    const bool DirBC_Interior = DirectSolver<LevelCacheType>::DirBC_Interior_;
+    const PolarGrid<DefaultMemorySpace>& grid = DirectSolver<LevelCacheType>::grid_;
+    const bool DirBC_Interior                 = DirectSolver<LevelCacheType>::DirBC_Interior_;
 
     assert(std::ssize(x) == grid.numberOfNodes());
     assert(grid.nr() >= 4);
@@ -108,4 +108,4 @@ void DirectSolverTake<LevelCacheType>::applySymmetryShift(HostVector<double> x)
     applySymmetryShiftOuterBoundary(x);
 
     Kokkos::fence();
-}
+}