SciCompMod · AbdelhadiKara · Apr 7, 2026 · Apr 7, 2026 · Apr 7, 2026 · Apr 8, 2026
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -24,7 +24,7 @@ if(NOT CMAKE_BUILD_TYPE)
 endif()
 
 # Set compiler flags
-set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -Wextra -pedantic -Wno-unused -Wno-psabi -Wfloat-conversion")
+set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fopenmp -Wall -Wextra -pedantic -Wno-unused -Wno-psabi -Wfloat-conversion")
 set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O2 -mtune=generic -Wno-psabi")
 
 # Set coverage compiler flags - must come before any targets are defined
@@ -44,6 +44,7 @@ if(GMGPOLAR_ENABLE_COVERAGE)
     endif()
 endif()
 
+find_package(Kokkos  4.4.1...<5.1 QUIET REQUIRED)
 include_directories(include)
 add_subdirectory(src)
 

diff --git a/include/PolarGrid/polargrid.h b/include/PolarGrid/polargrid.h
@@ -13,20 +13,22 @@
 #include <set>
 #include <stdexcept>
 #include <string>
-#include <vector>
 
 #include <omp.h>
 
 #include "../Definitions/equals.h"
+#include "../LinearAlgebra/Vector/vector.h"
 
 class PolarGrid
 {
 public:
     // Default constructor.
     explicit PolarGrid() = default;
 
-    // Constructor to initialize grid using vectors of radii and angles.
-    PolarGrid(const std::vector<double>& radii, const std::vector<double>& angles,
+    // Constructor to initialize grid using kokkos views of radii and angles.
+    PolarGrid(Vector<double> radii, Vector<double> angles, std::optional<double> splitting_radius = std::nullopt);
+    // Constructor to initialize grid using std::vectors of radii and angles.
+    PolarGrid(std::vector<double> radii, std::vector<double> angles,
               std::optional<double> splitting_radius = std::nullopt);
 
     // Constructor to initialize grid using parameters from GMGPolar.
@@ -76,19 +78,19 @@ class PolarGrid
     int nr_; // number of nodes in radial direction
     int ntheta_; // number of (unique) nodes in angular direction
     bool is_ntheta_PowerOfTwo_;
-    std::vector<double> radii_; // Vector of radial coordiantes
-    std::vector<double> angles_; // Vector of angular coordinates
+    AllocatableVector<double> radii_; // Vector of radial coordiantes
+    AllocatableVector<double> angles_; // Vector of angular coordinates
 
     // radial_spacings_ contains the distances between each consecutive radii division.
     // radial_spacings_ = [r_{1}-r_{0}, ..., r_{N}-r_{N-1}].
-    std::vector<double> radial_spacings_; // size(radial_spacings_) = nr() - 1
+    AllocatableVector<double> radial_spacings_; // size(radial_spacings_) = nr() - 1
 
     // angular_spacings_ contains the angles between each consecutive theta division.
     // Since we have a periodic boundary in theta direction,
     // we have to make sure the index wraps around correctly when accessing it.
     // Here theta_0 = 0.0 and theta_N = 2*pi refer to the same point.
     // angular_spacings_ = [theta_{1}-theta_{0}, ..., theta_{N}-theta_{N-1}].
-    std::vector<double> angular_spacings_; // size(angular_spacings_) = ntheta()
+    AllocatableVector<double> angular_spacings_; // size(angular_spacings_) = ntheta()
 
     // Circle/radial smoother division
     double smoother_splitting_radius_; // Radius at which the grid is split into circular and radial smoothing
@@ -105,7 +107,7 @@ class PolarGrid
      */
 
     // Check parameter validity
-    void checkParameters(const std::vector<double>& radii, const std::vector<double>& angles) const;
+    void checkParameters(Vector<double> radii, Vector<double> angles) const;
 
     // Initialize radial_spacings_, angular_spacings_
     void initializeDistances();
@@ -124,12 +126,12 @@ class PolarGrid
 
     // Refine the grid by dividing radial and angular divisions by 2.
     void refineGrid(const int divideBy2);
-    std::vector<double> divideVector(const std::vector<double>& vec, const int divideBy2) const;
+    Vector<double> divideVector(Vector<double> vec, const int divideBy2) const;
 
     // Help constrcut radii_ when an anisotropic radial division is requested
     // Implementation in src/PolarGrid/anisotropic_division.cpp
-    void RadialAnisotropicDivision(std::vector<double>& r_temp, double R0, double R, const int nr_exp,
-                                   double refinement_radius, const int anisotropic_factor) const;
+    Vector<double> RadialAnisotropicDivision(double R0, double R, const int nr_exp, double refinement_radius,
+                                             const int anisotropic_factor) const;
 };
 
 // ---------------------------------------------------- //

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
@@ -74,11 +74,10 @@ target_link_libraries(GMGPolarLib PUBLIC
 # Handle OpenMP
 find_package(OpenMP REQUIRED)
 if(OpenMP_CXX_FOUND)
-    target_link_libraries(GMGPolarLib PUBLIC OpenMP::OpenMP_CXX)
+    target_link_libraries(GMGPolarLib PUBLIC Kokkos::kokkos OpenMP::OpenMP_CXX)
+else()
+    target_link_libraries(GMGPolarLib PUBLIC Kokkos::kokkos)
 endif()
-find_package(Kokkos  4.4.1...<5.1 QUIET REQUIRED)
-target_link_libraries(GMGPolarLib PUBLIC Kokkos::kokkos)
-
 
 # Handle MUMPS configuration
 if(GMGPOLAR_USE_MUMPS)

diff --git a/src/PolarGrid/anisotropic_division.cpp b/src/PolarGrid/anisotropic_division.cpp
@@ -1,7 +1,7 @@
 #include "../../include/PolarGrid/polargrid.h"
 
-void PolarGrid::RadialAnisotropicDivision(std::vector<double>& r_temp, double R0, double R, const int nr_exp,
-                                          double refinement_radius, const int anisotropic_factor) const
+Vector<double> PolarGrid::RadialAnisotropicDivision(double R0, double R, const int nr_exp, double refinement_radius,
+                                                    const int anisotropic_factor) const
 {
     // Calculate the percentage of refinement_radius.
     const double percentage = (refinement_radius - R0) / (R - R0);
@@ -23,9 +23,9 @@ void PolarGrid::RadialAnisotropicDivision(std::vector<double>& r_temp, double R0
     if ((anisotropic_factor % 2) ==
         1) // odd number of elements on an open circular disk is desired because of coarsening
         n_elems_equi++;
-    double uniform_distance     = (R - R0) / n_elems_equi;
-    int nr                      = n_elems_equi + 1;
-    std::vector<double> r_temp2 = std::vector<double>(nr);
+    double uniform_distance = (R - R0) / n_elems_equi;
+    int nr                  = n_elems_equi + 1;
+    Vector<double> r_temp2("r_temp2", nr);
     for (int i = 0; i < nr - 1; i++)
         r_temp2[i] = R0 + i * uniform_distance;
     r_temp2[nr - 1] = R;
@@ -89,8 +89,7 @@ void PolarGrid::RadialAnisotropicDivision(std::vector<double>& r_temp, double R0
     // group all in r_tmp
     nr = n_elems_equi - n_elems_refined + r_set.size() + 1;
 
-    r_temp.resize(nr);
-
+    Vector<double> r_temp("r_temp", nr);
     for (int i = 0; i < se; i++)
         r_temp[i] = r_temp2[i];
     itr = r_set.begin();
@@ -100,4 +99,5 @@ void PolarGrid::RadialAnisotropicDivision(std::vector<double>& r_temp, double R0
     }
     for (int i = 0; i < n_elems_equi - ee + 1; i++)
         r_temp[se + r_set.size() + i] = r_temp2[ee + i];
+    return r_temp;
 }
diff --git a/src/PolarGrid/polargrid.cpp b/src/PolarGrid/polargrid.cpp
@@ -1,12 +1,11 @@
 #include "../../include/PolarGrid/polargrid.h"
-
+#include <Kokkos_StdAlgorithms.hpp>
 // ------------ //
 // Constructors //
 // ------------ //
 
 // Constructor to initialize grid using vectors of radii and angles.
-PolarGrid::PolarGrid(const std::vector<double>& radii, const std::vector<double>& angles,
-                     std::optional<double> splitting_radius)
+PolarGrid::PolarGrid(Vector<double> radii, Vector<double> angles, std::optional<double> splitting_radius)
     : nr_(radii.size())
     , ntheta_(angles.size() - 1)
     , is_ntheta_PowerOfTwo_((ntheta_ & (ntheta_ - 1)) == 0)
@@ -22,6 +21,32 @@ PolarGrid::PolarGrid(const std::vector<double>& radii, const std::vector<double>
     initializeLineSplitting(splitting_radius);
 }
 
+// Constructor to initialize grid using vectors of radii and angles.
+PolarGrid::PolarGrid(std::vector<double> radii, std::vector<double> angles, std::optional<double> splitting_radius)
+    : nr_(radii.size())
+    , ntheta_(angles.size() - 1)
+    , is_ntheta_PowerOfTwo_((ntheta_ & (ntheta_ - 1)) == 0)
+    , radii_("radii", nr_)
+    , angles_("angles", angles.size())
+
+{
+    // Copy from std vector to Kokkos view
+    for (std::size_t i(0); i < radii.size(); ++i) {
+        radii_(i) = radii[i];
+    }
+    for (std::size_t i(0); i < angles.size(); ++i) {
+        angles_(i) = angles[i];
+    }
+
+    // Check parameter validity
+    checkParameters(radii_, angles_);
+    // Store distances to adjacent neighboring nodes.
+    // Initializes radial_spacings_, angular_spacings_
+    initializeDistances();
+    // Initializes smoothers splitting radius for circle/radial indexing.
+    initializeLineSplitting(splitting_radius);
+}
+
 // Constructor to initialize grid using parameters from GMGPolar.
 PolarGrid::PolarGrid(double R0, double Rmax, const int nr_exp, const int ntheta_exp, double refinement_radius,
                      const int anisotropic_factor, const int divideBy2, std::optional<double> splitting_radius)
@@ -50,25 +75,25 @@ void PolarGrid::constructRadialDivisions(double R0, double R, const int nr_exp,
 {
     // r_temp contains the values before we refine one last time for extrapolation.
     // Therefore we first consider 2^(nr_exp-1) points.
-    std::vector<double> r_temp;
+    AllocatableVector<double> r_temp;
     if (anisotropic_factor == 0) {
         // nr = 2**(nr_exp-1) + 1
         int nr                  = (1 << (nr_exp - 1)) + 1;
         double uniform_distance = (R - R0) / (nr - 1);
         assert(uniform_distance > 0.0);
-        r_temp.resize(nr);
+        r_temp = Vector<double>("r_temp", nr);
         for (int i = 0; i < nr - 1; i++) {
             r_temp[i] = R0 + i * uniform_distance;
         }
         r_temp[nr - 1] = R;
     }
     else {
         // Implementation in src/PolarGrid/anisotropic_division.cpp
-        RadialAnisotropicDivision(r_temp, R0, R, nr_exp, refinement_radius, anisotropic_factor);
+        r_temp = RadialAnisotropicDivision(R0, R, nr_exp, refinement_radius, anisotropic_factor);
     }
     // Refine division in the middle for extrapolation
-    nr_ = 2 * r_temp.size() - 1;
-    radii_.resize(nr_);
+    nr_    = 2 * r_temp.size() - 1;
+    radii_ = AllocatableVector<double>("radii_", nr_);
     for (int i = 0; i < nr_; i++) {
         if (!(i % 2))
             radii_[i] = r_temp[i / 2];
@@ -94,7 +119,7 @@ void PolarGrid::constructAngularDivisions(const int ntheta_exp, const int nr)
     is_ntheta_PowerOfTwo_ = (ntheta_ & (ntheta_ - 1)) == 0;
     // Note that currently ntheta_ = 2^k which allows us to do some optimizations when indexing.
     double uniform_distance = 2 * M_PI / ntheta_;
-    angles_.resize(ntheta_ + 1);
+    Kokkos::resize(angles_, ntheta_ + 1);
     for (int i = 0; i < ntheta_; i++) {
         angles_[i] = i * uniform_distance;
     }
@@ -111,12 +136,12 @@ void PolarGrid::refineGrid(const int divideBy2)
     is_ntheta_PowerOfTwo_ = (ntheta_ & (ntheta_ - 1)) == 0;
 }
 
-std::vector<double> PolarGrid::divideVector(const std::vector<double>& vec, const int divideBy2) const
+Vector<double> PolarGrid::divideVector(Vector<double> vec, const int divideBy2) const
 {
     const int powerOfTwo = 1 << divideBy2;
     size_t vecSize       = vec.size();
     size_t resultSize    = vecSize + (vecSize - 1) * (powerOfTwo - 1);
-    std::vector<double> result(resultSize);
+    Vector<double> result("result", resultSize);
 
     for (size_t i = 0; i < vecSize - 1; ++i) {
         size_t baseIndex  = i * powerOfTwo;
@@ -126,15 +151,15 @@ std::vector<double> PolarGrid::divideVector(const std::vector<double>& vec, cons
             result[baseIndex + j]     = interpolated_value;
         }
     }
-    result[resultSize - 1] = vec.back(); // Add the last value of the original vector
+    result[resultSize - 1] = vec(vec.size() - 1); // Add the last value of the original vector
     return result;
 }
 
 void PolarGrid::initializeDistances()
 {
     // radial_spacings contains the distances between each consecutive radii division.
     // radial_spacings = [R_1-R0, ..., R_{N} - R_{N-1}].
-    radial_spacings_.resize(nr() - 1);
+    Kokkos::resize(radial_spacings_, nr() - 1);
     for (int i = 0; i < nr() - 1; i++) {
         radial_spacings_[i] = radius(i + 1) - radius(i);
     }
@@ -143,7 +168,7 @@ void PolarGrid::initializeDistances()
     // we have to make sure the index wraps around correctly when accessing it.
     // Here theta_0 = 0.0 and theta_N = 2*pi refer to the same point.
     // angular_spacings = [theta_{1}-theta_{0}, ..., theta_{N}-theta_{N-1}].
-    angular_spacings_.resize(ntheta());
+    Kokkos::resize(angular_spacings_, ntheta());
     for (int i = 0; i < ntheta(); i++) {
         angular_spacings_[i] = theta(i + 1) - theta(i);
     }
@@ -157,22 +182,25 @@ void PolarGrid::initializeDistances()
 void PolarGrid::initializeLineSplitting(std::optional<double> splitting_radius)
 {
     if (splitting_radius.has_value()) {
-        if (splitting_radius.value() < radii_.front()) {
+        if (splitting_radius.value() < radii_(0)) {
             number_smoother_circles_   = 0;
             length_smoother_radial_    = nr();
             smoother_splitting_radius_ = -1.0;
         }
         else {
-            auto it = std::lower_bound(radii_.begin(), radii_.end(), splitting_radius.value());
-            if (it != radii_.end()) {
-                number_smoother_circles_   = std::distance(radii_.begin(), it);
+            auto start = Kokkos::Experimental::begin(radii_);
+            auto end   = Kokkos::Experimental::end(radii_);
+            auto it    = std::lower_bound(start, end, splitting_radius.value());
+
+            if (it != end) {
+                number_smoother_circles_   = std::distance(start, it);
                 length_smoother_radial_    = nr() - number_smoother_circles_;
                 smoother_splitting_radius_ = splitting_radius.value();
             }
             else {
                 number_smoother_circles_   = nr();
                 length_smoother_radial_    = 0;
-                smoother_splitting_radius_ = radii_.back() + 1.0;
+                smoother_splitting_radius_ = radii_(radii_.size() - 1) + 1.0;
             }
         }
     }
@@ -215,8 +243,8 @@ PolarGrid coarseningGrid(const PolarGrid& fineGrid)
     const int coarse_nr     = (fineGrid.nr() + 1) / 2;
     const int coarse_ntheta = fineGrid.ntheta() / 2;
 
-    std::vector<double> coarse_r(coarse_nr);
-    std::vector<double> coarse_theta(coarse_ntheta + 1);
+    Vector<double> coarse_r("coarse_r", coarse_nr);
+    Vector<double> coarse_theta("coarse_theta", coarse_ntheta + 1);
 
     for (int i = 0; i < coarse_nr; i++) {
         coarse_r[i] = fineGrid.radius(2 * i);
@@ -239,41 +267,44 @@ PolarGrid coarseningGrid(const PolarGrid& fineGrid)
 // Check parameter validity //
 // ------------------------ //
 
-void PolarGrid::checkParameters(const std::vector<double>& radii, const std::vector<double>& angles) const
+void PolarGrid::checkParameters(Vector<double> radii, Vector<double> angles) const
 {
+    auto radii_start = Kokkos::Experimental::begin(radii);
+    auto radii_end   = Kokkos::Experimental::end(radii);
     if (radii.size() < 2) {
         throw std::invalid_argument("At least two radii are required.");
     }
 
-    if (!std::all_of(radii.begin(), radii.end(), [](double r) {
+    if (!std::all_of(radii_start, radii_end, [](double r) {
             return r > 0.0;
         })) {
         throw std::invalid_argument("All radii must be greater than zero.");
     }
 
-    if (std::adjacent_find(radii.begin(), radii.end(), std::greater_equal<double>()) != radii.end()) {
+    if (std::adjacent_find(radii_start, radii_end, std::greater_equal<double>()) != radii_end) {
         throw std::invalid_argument("Radii must be strictly increasing.");
     }
 
     if (angles.size() < 3) {
         throw std::invalid_argument("At least two angles are required.");
     }
-
-    if (!std::all_of(angles.begin(), angles.end(), [](double theta) {
+    auto angles_start = Kokkos::Experimental::begin(angles);
+    auto angles_end   = Kokkos::Experimental::end(angles);
+    if (!std::all_of(angles_start, angles_end, [](double theta) {
             return theta >= 0.0;
         })) {
         throw std::invalid_argument("All angles must be non-negative.");
     }
 
-    if (std::adjacent_find(angles.begin(), angles.end(), std::greater_equal<double>()) != angles.end()) {
+    if (std::adjacent_find(angles_start, angles_end, std::greater_equal<double>()) != angles_end) {
         throw std::invalid_argument("Angles must be strictly increasing.");
     }
 
-    if (!equals(angles.front(), 0.0)) {
+    if (!equals(angles(0), 0.0)) {
         throw std::invalid_argument("First angle must be 0.");
     }
 
-    if (!equals(angles.back(), 2 * M_PI)) {
+    if (!equals(angles(angles.size() - 1), 2 * M_PI)) {
         throw std::invalid_argument("Last angle must be 2*pi.");
     }
 

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
@@ -31,7 +31,7 @@ add_executable(gmgpolar_tests
 )
 
 # Set the compile features and link libraries
-target_compile_features(gmgpolar_tests PRIVATE cxx_std_17)
+target_compile_features(gmgpolar_tests PRIVATE cxx_std_20)
 target_link_libraries(gmgpolar_tests GMGPolarLib GTest::gtest_main)
 
 include(GoogleTest)

diff --git a/third-party/CMakeLists.txt b/third-party/CMakeLists.txt
@@ -4,7 +4,7 @@ include(FetchContent)
 FetchContent_Declare(
   googletest
   GIT_REPOSITORY https://github.com/google/googletest.git
-  GIT_TAG release-1.12.1
+  GIT_TAG v1.17.0
 )
 
 # For Windows: Prevent overriding the parent project's compiler/linker settings