[wip] feat: Add SYCL implementation of FTCS

CMakeLists.txt

@@ -99,4 +99,4 @@ endif()
 
 if(TUG_NAAICE_EXAMPLE)
   add_subdirectory(naaice)
-endif()
+endif()

examples/CMakeLists.txt

@@ -1,20 +1,29 @@
-add_executable(FTCS_1D_proto_example FTCS_1D_proto_example.cpp)
-add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp)
-add_executable(BTCS_1D_proto_example BTCS_1D_proto_example.cpp)
-add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp)
-add_executable(CRNI_2D_proto_example CRNI_2D_proto_example.cpp)
-add_executable(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp)
-add_executable(profiling_openmp profiling_openmp.cpp)
+# add_executable(FTCS_1D_proto_example FTCS_1D_proto_example.cpp)
+# add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp)
+# add_executable(BTCS_1D_proto_example BTCS_1D_proto_example.cpp)
+# add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp)
+# add_executable(CRNI_2D_proto_example CRNI_2D_proto_example.cpp)
+# add_executable(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp)
+# add_executable(profiling_openmp profiling_openmp.cpp)
 
-target_link_libraries(FTCS_1D_proto_example tug)
-target_link_libraries(FTCS_2D_proto_example tug)
-target_link_libraries(BTCS_1D_proto_example tug)
-target_link_libraries(BTCS_2D_proto_example tug)
-target_link_libraries(CRNI_2D_proto_example tug)
-target_link_libraries(reference-FTCS_2D_closed tug)
-target_link_libraries(profiling_openmp tug)
+# target_link_libraries(FTCS_1D_proto_example tug)
+# target_link_libraries(FTCS_2D_proto_example tug)
+# target_link_libraries(BTCS_1D_proto_example tug)
+# target_link_libraries(BTCS_2D_proto_example tug)
+# target_link_libraries(CRNI_2D_proto_example tug)
+# target_link_libraries(reference-FTCS_2D_closed tug)
+# target_link_libraries(profiling_openmp tug)
 
-add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp)
-add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp)
-target_link_libraries(FTCS_2D_proto_closed_mdl tug)
-target_link_libraries(FTCS_2D_proto_example_mdl tug)
+# add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp)
+# add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp)
+# target_link_libraries(FTCS_2D_proto_closed_mdl tug)
+# target_link_libraries(FTCS_2D_proto_example_mdl tug)
+
+find_package(IntelSYCL)
+
+if(IntelSYCL_FOUND)
+    add_executable(sycl_test FTCS_1D_sycl.cpp)
+    target_link_libraries(sycl_test PRIVATE tug)
+    target_compile_definitions(sycl_test PRIVATE TUG_ENABLE_SYCL)
+    add_sycl_to_target(TARGET sycl_test SOURCES FTCS_1D_sycl.cpp)
+endif()

examples/FTCS_1D_sycl.cpp

@@ -0,0 +1,93 @@
+#include <chrono>
+#include <tug/Grid.hpp>
+#include <tug/Simulation.hpp>
+
+#include <Eigen/Dense>
+
+using namespace tug;
+using namespace Eigen;
+
+MatrixXd runNormal(MatrixXd concentrations, const MatrixXd &alphas,
+                   Boundary<double> &boundaries, int iterations,
+                   double timestep) {
+  Grid64 grid(concentrations.cols());
+  grid.setConcentrations(concentrations);
+
+  grid.setAlpha(alphas);
+
+  Simulation simulation =
+      Simulation<double, tug::FTCS_APPROACH>(grid, boundaries);
+
+  simulation.setTimestep(timestep);
+  simulation.setIterations(iterations);
+  simulation.run();
+
+  return MatrixXd(grid.getConcentrations());
+}
+
+MatrixXd runWithSYCL(MatrixXd concentrations, const MatrixXd &alphas,
+                     Boundary<double> &boundaries, int iterations,
+                     double timestep) {
+  Grid64 grid(concentrations.cols());
+  grid.setConcentrations(concentrations);
+
+  grid.setAlpha(alphas);
+
+  Simulation simulation = Simulation<double, tug::FTCS_SYCL>(grid, boundaries);
+
+  simulation.setTimestep(timestep);
+  simulation.setIterations(iterations);
+  simulation.run();
+
+  return MatrixXd(grid.getConcentrations());
+}
+
+int main(int argc, char **argv) {
+  constexpr int cells = 1e7;
+  constexpr double timestep = 0.1;
+  constexpr int iterations = 500;
+
+  Grid64 grid(cells);
+
+  MatrixXd concentrations = MatrixXd::Constant(1, cells, 20);
+
+  MatrixXd alphas = MatrixXd::Constant(1, cells, 1);
+
+  // ******************
+  // **** BOUNDARY ****
+  // ******************
+
+  // create a boundary with constant values
+  Boundary bc = Boundary(grid);
+  bc.setBoundarySideConstant(BC_SIDE_LEFT, 1);
+  bc.setBoundarySideConstant(BC_SIDE_RIGHT, 1);
+
+  auto t_normal_begin = std::chrono::steady_clock::now();
+  auto normal = runNormal(concentrations, alphas, bc, iterations, timestep);
+  auto t_normal_end = std::chrono::steady_clock::now();
+
+  std::cout << "Normal time: "
+            << std::chrono::duration_cast<std::chrono::milliseconds>(
+                   t_normal_end - t_normal_begin)
+                   .count()
+            << "ms" << std::endl;
+
+  auto t_sycl_begin = std::chrono::steady_clock::now();
+  auto sycl = runWithSYCL(concentrations, alphas, bc, iterations, timestep);
+  auto t_sycl_end = std::chrono::steady_clock::now();
+
+  std::cout << "SYCL time: "
+            << std::chrono::duration_cast<std::chrono::milliseconds>(
+                   t_sycl_end - t_sycl_begin)
+                   .count()
+            << "ms" << std::endl;
+
+  // std::cout << "Normal: " << std::endl << normal << std::endl;
+  // std::cout << "SYCL: " << std::endl << sycl << std::endl;
+
+  auto diff = (normal - sycl).cwiseAbs().maxCoeff();
+
+  std::cout << "Max diff: " << diff << std::endl;
+
+  return 0;
+}

include/tug/Core/FTCS_Sycl.hpp

@@ -0,0 +1,139 @@
+#pragma once
+
+#include "../Boundary.hpp"
+#include "../Grid.hpp"
+#include "sycl/usm.hpp"
+
+#include <cstddef>
+#include <cstring>
+#include <sycl/sycl.hpp>
+
+namespace tug {
+
+template <typename T>
+static inline constexpr T boundApply(const T conc, const T conc_neighbor,
+                                     const T alpha, const T alpha_neighbor,
+                                     const T bcValue, const BC_TYPE bcType) {
+  const T alpha_intercell = calcAlphaIntercell(alpha, alpha_neighbor);
+  switch (bcType) {
+  case BC_TYPE::BC_TYPE_CLOSED: {
+    const T conc_change = conc_neighbor - conc;
+
+    return alpha_intercell * conc_change;
+  }
+  case BC_TYPE::BC_TYPE_CONSTANT: {
+    const T inner_term = alpha_intercell * conc_neighbor;
+    const T mid_term = (alpha_intercell + 2 * alpha) * conc;
+    const T outer_term = 2 * alpha * bcValue;
+
+    return inner_term - mid_term + outer_term;
+  }
+  default:;
+  }
+}
+
+template <typename T>
+static inline constexpr T explicitChange(const T conc, const T conc_n1,
+                                         const T conc_n2, const T alpha,
+                                         const T alpha_n1, const T alpha_n2) {
+  const T alpha_intercell_n1 = calcAlphaIntercell(alpha, alpha_n1);
+  const T alpha_intercell_n2 = calcAlphaIntercell(alpha, alpha_n2);
+
+  const T conc_term_n1 = alpha_intercell_n1 * conc_n1;
+  const T conc_term_n2 = alpha_intercell_n2 * conc_n2;
+
+  return conc_term_n1 + conc_term_n2 -
+         (alpha_intercell_n1 + alpha_intercell_n2) * conc;
+}
+
+template <typename T>
+static inline void
+FTCS_SYCL_1D(sycl::queue &q, T *concentrations_inout, const T *alphas_in,
+             const tug::BoundaryElement<T> &leftbound,
+             const tug::BoundaryElement<T> &rightbound,
+             const std::size_t colMax, const T dtOverDeltaColSquare,
+             const std::size_t nIter) {
+
+  const BC_TYPE bcType_left = leftbound.getType();
+  const BC_TYPE bcType_right = rightbound.getType();
+
+  const T bcValue_left = leftbound.getValue();
+  const T bcValue_right = rightbound.getValue();
+
+  {
+    T *concentrations_device = sycl::malloc_device<T>(colMax, q);
+    T *concentrations_t1_device = sycl::malloc_device<T>(colMax, q);
+    T *alpha_device = sycl::malloc_device<T>(colMax, q);
+
+    q.memcpy(concentrations_device, concentrations_inout, colMax * sizeof(T));
+    q.memcpy(alpha_device, alphas_in, colMax * sizeof(T));
+
+    q.wait();
+
+    for (std::size_t i = 0; i < nIter; i++) {
+      auto e_boundleft = q.submit([&](sycl::handler &h) {
+        h.single_task([=]() {
+          concentrations_t1_device[0] =
+              concentrations_device[0] +
+              dtOverDeltaColSquare *
+                  boundApply(concentrations_device[0], concentrations_device[1],
+                             alpha_device[0], alpha_device[1], bcValue_left,
+                             bcType_left);
+        });
+      });
+
+      auto e_boundright = q.submit([&](sycl::handler &h) {
+        h.single_task([=]() {
+          concentrations_t1_device[colMax - 1] =
+              concentrations_device[colMax - 1] +
+              dtOverDeltaColSquare *
+                  boundApply(concentrations_device[colMax - 1],
+                             concentrations_device[colMax - 2],
+                             alpha_device[colMax - 1], alpha_device[colMax - 2],
+                             bcValue_right, bcType_right);
+        });
+      });
+
+      auto e_mid = q.submit([&](sycl::handler &h) {
+        h.parallel_for(sycl::range<1>(colMax - 2), [=](sycl::id<1> idx) {
+          const std::size_t i = idx[0] + 1;
+          concentrations_t1_device[i] =
+              concentrations_device[i] +
+              dtOverDeltaColSquare *
+                  explicitChange(concentrations_device[i],
+                                 concentrations_device[i + 1],
+                                 concentrations_device[i - 1], alpha_device[i],
+                                 alpha_device[i + 1], alpha_device[i - 1]);
+        });
+      });
+
+      q.memcpy(concentrations_device, concentrations_t1_device,
+               colMax * sizeof(T), {e_boundleft, e_boundright, e_mid})
+          .wait();
+    }
+
+    q.memcpy(concentrations_inout, concentrations_device, colMax * sizeof(T))
+        .wait();
+    sycl::free(concentrations_device, q);
+    sycl::free(concentrations_t1_device, q);
+    sycl::free(alpha_device, q);
+  }
+}
+
+template <typename T>
+void FTCS_Sycl_impl(sycl::queue &q, tug::Grid<T> &grid,
+                    const tug::Boundary<T> &boundaries, const T dt,
+                    const std::size_t nIter) {
+  const std::size_t colMax = grid.getCol();
+  const T deltaCol = grid.getDeltaCol();
+  const T deltaColSquare = deltaCol * deltaCol;
+  const T dtOverDeltaColSquare = dt / deltaColSquare;
+
+  if (grid.getDim() == 1) {
+    FTCS_SYCL_1D(q, grid.getConcentrations().data(), grid.getAlpha().data(),
+                 boundaries.getBoundaryElement(BC_SIDE_LEFT, 0),
+                 boundaries.getBoundaryElement(BC_SIDE_RIGHT, 0), colMax,
+                 dtOverDeltaColSquare, nIter);
+  }
+}
+} // namespace tug

include/tug/Grid.hpp

@@ -227,15 +227,7 @@ public:
    *
    * @return A matrix holding the alpha coefficients in x-direction.
    */
-  const auto &getAlphaX() const {
-
-    if (dim != 2) {
-      throw std::invalid_argument(
-          "Grid is not two dimensional, you should probably use getAlpha()!");
-    }
-
-    return this->alphaX;
-  }
+  const auto &getAlphaX() const { return this->alphaX; }
 
   /**
    * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid.

include/tug/Simulation.hpp

@@ -23,6 +23,11 @@
 #include "Core/FTCS.hpp"
 #include "Core/TugUtils.hpp"
 
+#ifdef TUG_ENABLE_SYCL
+#include "Core/FTCS_Sycl.hpp"
+#include <sycl/sycl.hpp>
+#endif
+
 #ifdef _OPENMP
 #include <omp.h>
 #else
@@ -39,6 +44,10 @@ enum APPROACH {
   FTCS_APPROACH,          /*!< Forward Time-Centered Space */
   BTCS_APPROACH,          /*!< Backward Time-Centered Space */
   CRANK_NICOLSON_APPROACH /*!< Crank-Nicolson method */
+#ifdef TUG_ENABLE_SYCL
+  ,
+  FTCS_SYCL /*!< Forward Time-Centered Space with SYCL */
+#endif
 };
 
 /**
@@ -108,7 +117,13 @@ public:
    * @param bc Valid boundary condition object
    * @param approach Approach to solving the problem. Either FTCS or BTCS.
    */
-  Simulation(Grid<T> &_grid, Boundary<T> &_bc) : grid(_grid), bc(_bc){};
+  Simulation(Grid<T> &_grid, Boundary<T> &_bc) : grid(_grid), bc(_bc) {
+#ifdef TUG_ENABLE_SYCL
+    if constexpr (approach == FTCS_SYCL) {
+      q = std::make_unique<sycl::queue>(sycl::default_selector_v);
+    }
+#endif
+  };
 
   /**
    * @brief Set the option to output the results to a CSV file. Off by default.
@@ -182,8 +197,7 @@ public:
       throw_invalid_argument("Timestep has to be greater than zero.");
     }
 
-    if constexpr (approach == FTCS_APPROACH ||
-                  approach == CRANK_NICOLSON_APPROACH) {
+    if constexpr (approach != BTCS_APPROACH) {
       T cfl;
       if (grid.getDim() == 1) {
 
@@ -412,7 +426,14 @@ public:
         // }
       }
 
-    } else if constexpr (approach == BTCS_APPROACH) { // BTCS case
+    }
+#ifdef TUG_ENABLE_SYCL
+    else if constexpr (approach == FTCS_SYCL) {
+      FTCS_Sycl_impl(*this->q.get(), this->grid, this->bc, this->timestep,
+                     this->iterations * innerIterations);
+    }
+#endif
+    else if constexpr (approach == BTCS_APPROACH) { // BTCS case
 
       if constexpr (solver == EIGEN_LU_SOLVER) {
         for (int i = 0; i < iterations; i++) {
@@ -498,7 +519,16 @@ private:
   Grid<T> &grid;
   Boundary<T> &bc;
 
-  const std::vector<std::string> approach_names = {"FTCS", "BTCS", "CRNI"};
+  const std::vector<std::string> approach_names = {"FTCS", "BTCS", "CRNI"
+#ifdef TUG_ENABLE_SYCL
+                                                   ,
+                                                   "FTCS_SYCL"
+#endif
+  };
+
+#ifdef TUG_ENABLE_SYCL
+  std::unique_ptr<sycl::queue> q;
+#endif
 };
 } // namespace tug
 #endif // SIMULATION_H_