Merge pull request #97 from Exawind/refactor_for_cuda_run

Fixed runtime errors when running unit tests in a CUDA build

src/CMakeLists.txt

@@ -2,15 +2,23 @@
 set(Kokkos_DIR "$ENV{Kokkos_ROOT}" CACHE STRING "Kokkos root directory")
 find_package(KokkosKernels REQUIRED)
 
+if(CMAKE_CXX_COMPILER_ID MATCHES GNU)
+  if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "9.0")
+    set(FS_LIB stdc++fs)
+  endif()
+endif()
+
 target_sources(${oturb_exe_name} PRIVATE main.cpp)
 set_target_properties(${oturb_exe_name} PROPERTIES LINKER_LANGUAGE CXX)
 target_link_libraries(${oturb_exe_name} PRIVATE
     Kokkos::kokkoskernels
+    ${FS_LIB}
 )
 
 set_target_properties(${oturb_lib_name} PROPERTIES LINKER_LANGUAGE CXX)
 target_link_libraries(${oturb_lib_name} PRIVATE
     Kokkos::kokkoskernels
+    ${FS_LIB}
 )
 
 target_include_directories(${oturb_exe_name} PRIVATE ${PROJECT_BINARY_DIR})

src/gebt_poc/element.cpp

@@ -1,5 +1,7 @@
 #include "src/gebt_poc/element.h"
 
+#include <KokkosBlas.hpp>
+
 #include "src/utilities/log.h"
 
 namespace openturbine::gebt_poc {
@@ -38,4 +40,14 @@ double CalculateJacobian(
     return jacobian;
 }
 
+double CalculateJacobian(
+    Kokkos::View<const double* [3]> nodes, Kokkos::View<const double*> shape_derivatives
+) {
+    auto v0 = KokkosBlas::dot(shape_derivatives, Kokkos::subview(nodes, Kokkos::ALL, 0));
+    auto v1 = KokkosBlas::dot(shape_derivatives, Kokkos::subview(nodes, Kokkos::ALL, 1));
+    auto v2 = KokkosBlas::dot(shape_derivatives, Kokkos::subview(nodes, Kokkos::ALL, 2));
+
+    return std::sqrt(v0 * v0 + v1 * v1 + v2 * v2);
+}
+
 }  // namespace openturbine::gebt_poc

src/gebt_poc/element.h

@@ -9,7 +9,7 @@ namespace openturbine::gebt_poc {
 // Physically, the Jacobian represents the following (In 1-D):
 // Length of element in the global csys / Length of element in the natural csys
 double CalculateJacobian(
-    const std::vector<Point>& nodes, const std::vector<double>& shape_derivatives
+    Kokkos::View<const double* [3]> nodes, Kokkos::View<const double*> shape_derivatives
 );
 
 }  // namespace openturbine::gebt_poc

src/gebt_poc/mesh.h

@@ -10,7 +10,7 @@ class Mesh {
     int GetNumberOfNodes() const { return number_of_nodes_; }
 
     KOKKOS_FUNCTION
-    Kokkos::View<const int*> GetNodesForElement(int element_id) const {
+    auto GetNodesForElement(int element_id) const {
         return Kokkos::subview(connectivity_, element_id, Kokkos::ALL);
     }
 
@@ -55,8 +55,10 @@ class Mesh {
     void CheckNodeUniquenessConsistency() {
         for (std::size_t element = 0; element < connectivity_.extent(0); ++element) {
             auto node_list = GetNodesForElement(element);
-            auto host_node_list = Kokkos::create_mirror(node_list);
-            Kokkos::deep_copy(host_node_list, node_list);
+            auto tmp_node_list = Kokkos::View<double*>("tmp", node_list.extent(0));
+            auto host_node_list = Kokkos::create_mirror(tmp_node_list);
+            Kokkos::deep_copy(tmp_node_list, node_list);
+            Kokkos::deep_copy(host_node_list, tmp_node_list);
             auto node_list_vector = std::vector<int>(host_node_list.extent(0));
             for (std::size_t node = 0; node < node_list_vector.size(); ++node) {
                 node_list_vector[node] = host_node_list(node);

src/gebt_poc/solver.cpp

@@ -29,7 +29,7 @@ void Interpolate(
                 value += interpolation_function(j) *
                          nodal_values(j * kNumberOfLieAlgebraComponents + i) / jacobian;
             },
-            interpolated_values(i)
+            Kokkos::Sum(Kokkos::subview(interpolated_values, i))
         );
     }
 }
@@ -39,8 +39,9 @@ void CalculateCurvature(
     Kokkos::View<double*> curvature
 ) {
     // curvature = B * q_prime
-    auto b_matrix = gen_alpha_solver::BMatrixForQuaternions(
-        gen_alpha_solver::Quaternion(gen_coords(3), gen_coords(4), gen_coords(5), gen_coords(6))
+    auto b_matrix = Kokkos::View<double[3][4]>("b_matrix");
+    gen_alpha_solver::BMatrixForQuaternions(
+        b_matrix, Kokkos::subview(gen_coords, Kokkos::make_pair(3, 7))
     );
     auto q_prime = Kokkos::subview(gen_coords_derivative, Kokkos::make_pair(3, 7));
     KokkosBlas::gemv("N", 2., b_matrix, q_prime, 0., curvature);
@@ -150,12 +151,12 @@ void CalculateStaticResidual(
     const auto order = n_nodes - 1;
     const auto n_quad_pts = quadrature.GetNumberOfQuadraturePoints();
 
-    std::vector<Point> nodes;
-    for (size_t i = 0; i < n_nodes; ++i) {
-        nodes.emplace_back(
-            position_vectors(i * kNumberOfLieAlgebraComponents),
-            position_vectors(i * kNumberOfLieAlgebraComponents + 1),
-            position_vectors(i * kNumberOfLieAlgebraComponents + 2)
+    auto nodes = Kokkos::View<double* [3]>("nodes", n_nodes);
+    for (std::size_t i = 0; i < n_nodes; ++i) {
+        auto index = i * kNumberOfLieAlgebraComponents;
+        Kokkos::deep_copy(
+            Kokkos::subview(nodes, i, Kokkos::ALL),
+            Kokkos::subview(position_vectors, Kokkos::make_pair(index, index + 3))
         );
     }
 
@@ -184,7 +185,7 @@ void CalculateStaticResidual(
             auto shape_function_derivative =
                 gen_alpha_solver::create_vector(LagrangePolynomialDerivative(order, q_pt));
 
-            auto jacobian = CalculateJacobian(nodes, LagrangePolynomialDerivative(order, q_pt));
+            auto jacobian = CalculateJacobian(nodes, shape_function_derivative);
             Interpolate(gen_coords, shape_function, 1., gen_coords_qp);
             Interpolate(gen_coords, shape_function_derivative, jacobian, gen_coords_derivatives_qp);
             Interpolate(position_vectors, shape_function, 1., position_vector_qp);
@@ -364,12 +365,12 @@ void CalculateStaticIterationMatrix(
     const auto order = n_nodes - 1;
     const auto n_quad_pts = quadrature.GetNumberOfQuadraturePoints();
 
-    std::vector<Point> nodes;
-    for (size_t i = 0; i < n_nodes; ++i) {
-        nodes.emplace_back(
-            position_vectors(i * kNumberOfLieAlgebraComponents),
-            position_vectors(i * kNumberOfLieAlgebraComponents + 1),
-            position_vectors(i * kNumberOfLieAlgebraComponents + 2)
+    auto nodes = Kokkos::View<double* [3]>("nodes", n_nodes);
+    for (std::size_t i = 0; i < n_nodes; ++i) {
+        auto index = i * kNumberOfLieAlgebraComponents;
+        Kokkos::deep_copy(
+            Kokkos::subview(nodes, i, Kokkos::ALL),
+            Kokkos::subview(position_vectors, Kokkos::make_pair(index, index + 3))
         );
     }
 
@@ -405,7 +406,7 @@ void CalculateStaticIterationMatrix(
                 auto shape_function_derivative =
                     gen_alpha_solver::create_vector(LagrangePolynomialDerivative(order, q_pt));
 
-                auto jacobian = CalculateJacobian(nodes, LagrangePolynomialDerivative(order, q_pt));
+                auto jacobian = CalculateJacobian(nodes, shape_function_derivative);
                 Interpolate(gen_coords, shape_function, 1., gen_coords_qp);
                 Interpolate(
                     gen_coords, shape_function_derivative, jacobian, gen_coords_derivatives_qp
@@ -539,11 +540,14 @@ void ConstraintsGradientMatrix(
         constraints_gradient_matrix, Kokkos::make_pair(3, 6), Kokkos::make_pair(3, 6)
     );
 
-    auto hostB11 = Kokkos::create_mirror_view(B11);
-    hostB11(0, 0) = 1.;
-    hostB11(1, 1) = 1.;
-    hostB11(2, 2) = 1.;
-    Kokkos::deep_copy(B11, hostB11);
+    Kokkos::parallel_for(
+        1,
+        KOKKOS_LAMBDA(std::size_t) {
+            B11(0, 0) = 1.;
+            B11(1, 1) = 1.;
+            B11(2, 2) = 1.;
+        }
+    );
 
     KokkosBlas::scal(B21, -1., rotation_matrix_0);
     KokkosBlas::gemm("N", "N", -1., rotation_matrix_0, position_cross_prod_matrix, 0., B22);

src/gen_alpha_poc/quaternion.h

@@ -249,21 +249,25 @@ RotationMatrix quaternion_to_rotation_matrix(const Quaternion& quaternion) {
 /// Converts a 4x1 unit quaternion to a 3x3 rotation matrix and returns the result
 inline Kokkos::View<double**> EulerParameterToRotationMatrix(const Kokkos::View<double*> euler_param
 ) {
-    auto c0 = euler_param(0);
-    auto c = Kokkos::View<double*>("c", 3);
+    auto c = Kokkos::View<double[3]>("c");
     Kokkos::parallel_for(
         3, KOKKOS_LAMBDA(const size_t i) { c(i) = euler_param(i + 1); }
     );
     auto identity_matrix = gen_alpha_solver::create_identity_matrix(3);
     auto tilde_c = gen_alpha_solver::create_cross_product_matrix(c);
     auto tilde_c_tilde_c = gen_alpha_solver::multiply_matrix_with_matrix(tilde_c, tilde_c);
 
-    auto rotation_matrix = Kokkos::View<double**>("rotation_matrix", 3, 3);
+    auto rotation_matrix = Kokkos::View<double[3][3]>("rotation_matrix");
     Kokkos::parallel_for(
-        Kokkos::MDRangePolicy<Kokkos::DefaultExecutionSpace, Kokkos::Rank<2>>({0, 0}, {3, 3}),
-        KOKKOS_LAMBDA(const size_t i, const size_t j) {
-            rotation_matrix(i, j) =
-                identity_matrix(i, j) + 2 * c0 * tilde_c(i, j) + 2 * tilde_c_tilde_c(i, j);
+        1,
+        KOKKOS_LAMBDA(std::size_t) {
+            for (std::size_t i = 0; i < 3; ++i) {
+                for (std::size_t j = 0; j < 3; ++j) {
+                    rotation_matrix(i, j) = identity_matrix(i, j) +
+                                            2 * euler_param(0) * tilde_c(i, j) +
+                                            2 * tilde_c_tilde_c(i, j);
+                }
+            }
         }
     );
     return rotation_matrix;
@@ -308,32 +312,26 @@ Quaternion rotation_matrix_to_quaternion(const RotationMatrix& rotation_matrix)
 }
 
 /// Returns the B derivative matrix given for Euler parameters, i.e. unit quaternions
-inline Kokkos::View<double**> BMatrixForQuaternions(const Quaternion& quaternion) {
-    auto q0 = quaternion.GetScalarComponent();
-    auto q1 = quaternion.GetXComponent();
-    auto q2 = quaternion.GetYComponent();
-    auto q3 = quaternion.GetZComponent();
-
-    Kokkos::View<double**> bmatrix("bmatrix", 3, 4);
+inline void BMatrixForQuaternions(
+    Kokkos::View<double[3][4]> bmatrix, Kokkos::View<const double[4]> quaternion
+) {
     auto populate_bmatrix = KOKKOS_LAMBDA(size_t) {
-        bmatrix(0, 0) = -q1;
-        bmatrix(0, 1) = q0;
-        bmatrix(0, 2) = -q3;
-        bmatrix(0, 3) = q2;
-
-        bmatrix(1, 0) = -q2;
-        bmatrix(1, 1) = q3;
-        bmatrix(1, 2) = q0;
-        bmatrix(1, 3) = -q1;
-
-        bmatrix(2, 0) = -q3;
-        bmatrix(2, 1) = -q2;
-        bmatrix(2, 2) = q1;
-        bmatrix(2, 3) = q0;
+        bmatrix(0, 0) = -quaternion(1);
+        bmatrix(0, 1) = quaternion(0);
+        bmatrix(0, 2) = -quaternion(3);
+        bmatrix(0, 3) = quaternion(2);
+
+        bmatrix(1, 0) = -quaternion(2);
+        bmatrix(1, 1) = quaternion(3);
+        bmatrix(1, 2) = quaternion(0);
+        bmatrix(1, 3) = -quaternion(1);
+
+        bmatrix(2, 0) = -quaternion(3);
+        bmatrix(2, 1) = -quaternion(2);
+        bmatrix(2, 2) = quaternion(1);
+        bmatrix(2, 3) = quaternion(0);
     };
     Kokkos::parallel_for(1, populate_bmatrix);
-
-    return bmatrix;
 }
 
 }  // namespace openturbine::gen_alpha_solver

src/gen_alpha_poc/utilities.cpp

@@ -257,21 +257,4 @@ Kokkos::View<double**> add_matrix_with_matrix(
     return add_matrix_with_matrix(matrix_a_copy, matrix_b_copy);
 }
 
-Kokkos::View<double[1]> dot_product(
-    const Kokkos::View<double*> vector_a, const Kokkos::View<double*> vector_b
-) {
-    if (vector_a.extent(0) != vector_b.extent(0)) {
-        throw std::invalid_argument("The dimensions of the vectors must be equal to each other");
-    }
-
-    auto result = Kokkos::View<double[1]>("result", 1);
-    auto entries = Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, vector_a.extent(0));
-    auto dot_product = KOKKOS_LAMBDA(size_t index) {
-        result(0) += vector_a(index) * vector_b(index);
-    };
-    Kokkos::parallel_for(entries, dot_product);
-
-    return result;
-}
-
 }  // namespace openturbine::gen_alpha_solver

src/gen_alpha_poc/utilities.h

@@ -82,7 +82,4 @@ add_matrix_with_matrix(const Kokkos::View<double**>, const Kokkos::View<double**
 Kokkos::View<double**>
 add_matrix_with_matrix(const Kokkos::View<const double**>, const Kokkos::View<const double**>);
 
-/// Multiplies an 1 x n vector with an n x 1 vector and returns a double
-Kokkos::View<double[1]> dot_product(const Kokkos::View<double*>, const Kokkos::View<double*>);
-
 }  // namespace openturbine::gen_alpha_solver

tests/unit_tests/CMakeLists.txt

@@ -18,9 +18,15 @@ target_compile_options(
 target_include_directories(${oturb_unit_test_exe_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
 target_include_directories(${oturb_unit_test_exe_name} PRIVATE ${PROJECT_BINARY_DIR})
 
+if(CMAKE_CXX_COMPILER_ID MATCHES GNU)
+  if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "9.0")
+    set(FS_LIB stdc++fs)
+  endif()
+endif()
+
 # Link our unit test executable with GoogleTest
 find_package(GTest REQUIRED)
-target_link_libraries(${oturb_unit_test_exe_name} PRIVATE GTest::gtest GTest::gtest_main)
+target_link_libraries(${oturb_unit_test_exe_name} PRIVATE ${FS_LIB} GTest::gtest GTest::gtest_main)
 
 # Link to OpenTurbine test targets
 target_link_libraries(${oturb_unit_test_exe_name} PRIVATE ${oturb_lib_name})

tests/unit_tests/gebt_poc/test_element.cpp

@@ -1,17 +1,33 @@
+#include <initializer_list>
+
 #include <gtest/gtest.h>
 
 #include "src/gebt_poc/element.h"
 #include "tests/unit_tests/gen_alpha_poc/test_utilities.h"
 
 namespace openturbine::gebt_poc::tests {
 
+auto MakePointVector(std::initializer_list<Point> point_list) {
+    auto points = Kokkos::View<double* [3]>("points", point_list.size());
+    auto host_points = Kokkos::create_mirror(points);
+    std::size_t i = 0;
+    for (auto& point : point_list) {
+        host_points(i, 0) = point.GetXComponent();
+        host_points(i, 1) = point.GetYComponent();
+        host_points(i, 2) = point.GetZComponent();
+        ++i;
+    }
+    Kokkos::deep_copy(points, host_points);
+    return points;
+}
+
 TEST(ElementTest, JacobianFor1DLinearElement) {
-    auto shape_derivatives = LagrangePolynomialDerivative(1, -1.);
+    auto shape_derivatives = gen_alpha_solver::create_vector(LagrangePolynomialDerivative(1, -1.));
 
     // Nodes are located at GLL points i.e. (-1, 0., 0.) and (1., 0., 0.)
     auto node_1 = Point(-1., 0., 0.);
     auto node_2 = Point(1., 0., 0.);
-    auto nodes_1 = std::vector<Point>{node_1, node_2};
+    auto nodes_1 = MakePointVector({node_1, node_2});
     auto jacobian_1 = CalculateJacobian(nodes_1, shape_derivatives);
 
     // ** For 1-D elements, the Jacobian is always length of the element / 2 **
@@ -20,34 +36,38 @@ TEST(ElementTest, JacobianFor1DLinearElement) {
 
     auto node_3 = Point(-2., 0., 0.);
     auto node_4 = Point(2., 0., 0.);
-    auto nodes_2 = std::vector<Point>{node_3, node_4};
+    auto nodes_2 = MakePointVector({node_3, node_4});
     auto jacobian_2 = CalculateJacobian(nodes_2, shape_derivatives);
 
     // Length of the element is 4, so Jacobian is 2
     EXPECT_DOUBLE_EQ(jacobian_2, 2.);
 
     auto node_5 = Point(0., 0., 0.);
     auto node_6 = Point(1., 0., 0.);
-    auto nodes_3 = std::vector<Point>{node_5, node_6};
+    auto nodes_3 = MakePointVector({node_5, node_6});
     auto jacobian_3 = CalculateJacobian(nodes_3, shape_derivatives);
 
     // Length of the element is 1, so Jacobian is 0.5
     EXPECT_DOUBLE_EQ(jacobian_3, 0.5);
 }
 
 TEST(ElementTest, JacobianFor1DFourthOrderElement) {
-    auto shape_derivatives_1 = LagrangePolynomialDerivative(4, -0.9061798459386640);
-    auto shape_derivatives_2 = LagrangePolynomialDerivative(4, -0.5384693101056831);
-    auto shape_derivatives_3 = LagrangePolynomialDerivative(4, 0.);
-    auto shape_derivatives_4 = LagrangePolynomialDerivative(4, 0.5384693101056831);
-    auto shape_derivatives_5 = LagrangePolynomialDerivative(4, 0.9061798459386640);
+    auto shape_derivatives_1 =
+        gen_alpha_solver::create_vector(LagrangePolynomialDerivative(4, -0.9061798459386640));
+    auto shape_derivatives_2 =
+        gen_alpha_solver::create_vector(LagrangePolynomialDerivative(4, -0.5384693101056831));
+    auto shape_derivatives_3 = gen_alpha_solver::create_vector(LagrangePolynomialDerivative(4, 0.));
+    auto shape_derivatives_4 =
+        gen_alpha_solver::create_vector(LagrangePolynomialDerivative(4, 0.5384693101056831));
+    auto shape_derivatives_5 =
+        gen_alpha_solver::create_vector(LagrangePolynomialDerivative(4, 0.9061798459386640));
 
     auto node_1 = Point(0.0, 0.0, 0.0);
     auto node_2 = Point(0.16237631096713473, 0.17578464768961147, 0.1481911137890286);
     auto node_3 = Point(0.25, 1., 1.1875);
     auto node_4 = Point(-0.30523345382427747, 2.4670724951675314, 2.953849702537502);
     auto node_5 = Point(-1., 3.5, 4.);
-    auto nodes = std::vector<Point>{node_1, node_2, node_3, node_4, node_5};
+    auto nodes = MakePointVector({node_1, node_2, node_3, node_4, node_5});
 
     auto jacobian_1 = CalculateJacobian(nodes, shape_derivatives_1);
     auto jacobian_2 = CalculateJacobian(nodes, shape_derivatives_2);