diff --git a/CMakeLists.txt b/CMakeLists.txt index 78e45c9..1f04574 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -99,4 +99,4 @@ endif() if(TUG_NAAICE_EXAMPLE) add_subdirectory(naaice) -endif() +endif() \ No newline at end of file diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 0a3173d..4794a49 100644 --- a/examples/CMakeLists.txt +++ b/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() \ No newline at end of file diff --git a/examples/FTCS_1D_sycl.cpp b/examples/FTCS_1D_sycl.cpp new file mode 100644 index 0000000..eb688cb --- /dev/null +++ b/examples/FTCS_1D_sycl.cpp @@ -0,0 +1,93 @@ +#include +#include +#include + +#include + +using namespace tug; +using namespace Eigen; + +MatrixXd runNormal(MatrixXd concentrations, const MatrixXd &alphas, + Boundary &boundaries, int iterations, + double timestep) { + Grid64 grid(concentrations.cols()); + grid.setConcentrations(concentrations); + + grid.setAlpha(alphas); + + Simulation simulation = + Simulation(grid, boundaries); + + simulation.setTimestep(timestep); + simulation.setIterations(iterations); + simulation.run(); + + return MatrixXd(grid.getConcentrations()); +} + +MatrixXd runWithSYCL(MatrixXd concentrations, const MatrixXd &alphas, + Boundary &boundaries, int iterations, + double timestep) { + Grid64 grid(concentrations.cols()); + grid.setConcentrations(concentrations); + + grid.setAlpha(alphas); + + Simulation simulation = Simulation(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( + 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( + 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; +} \ No newline at end of file diff --git a/include/tug/Core/FTCS_Sycl.hpp b/include/tug/Core/FTCS_Sycl.hpp new file mode 100644 index 0000000..28479c5 --- /dev/null +++ b/include/tug/Core/FTCS_Sycl.hpp @@ -0,0 +1,139 @@ +#pragma once + +#include "../Boundary.hpp" +#include "../Grid.hpp" +#include "sycl/usm.hpp" + +#include +#include +#include + +namespace tug { + +template +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 +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 +static inline void +FTCS_SYCL_1D(sycl::queue &q, T *concentrations_inout, const T *alphas_in, + const tug::BoundaryElement &leftbound, + const tug::BoundaryElement &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(colMax, q); + T *concentrations_t1_device = sycl::malloc_device(colMax, q); + T *alpha_device = sycl::malloc_device(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 +void FTCS_Sycl_impl(sycl::queue &q, tug::Grid &grid, + const tug::Boundary &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 \ No newline at end of file diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index b8e471a..2e2c0de 100644 --- a/include/tug/Grid.hpp +++ b/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. diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index 7f8632d..549d7da 100644 --- a/include/tug/Simulation.hpp +++ b/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 +#endif + #ifdef _OPENMP #include #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 &_grid, Boundary &_bc) : grid(_grid), bc(_bc){}; + Simulation(Grid &_grid, Boundary &_bc) : grid(_grid), bc(_bc) { +#ifdef TUG_ENABLE_SYCL + if constexpr (approach == FTCS_SYCL) { + q = std::make_unique(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 &grid; Boundary &bc; - const std::vector approach_names = {"FTCS", "BTCS", "CRNI"}; + const std::vector approach_names = {"FTCS", "BTCS", "CRNI" +#ifdef TUG_ENABLE_SYCL + , + "FTCS_SYCL" +#endif + }; + +#ifdef TUG_ENABLE_SYCL + std::unique_ptr q; +#endif }; } // namespace tug #endif // SIMULATION_H_