From eadbd325ed69abf148351e52c0a270ab32597169 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Sat, 20 Dec 2025 00:02:45 +0100 Subject: Working example of poiseulle on gpu --- examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp | 113 +++++++++++++++++-------- 1 file changed, 78 insertions(+), 35 deletions(-) diff --git a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp index b31f8c5..2c50dd4 100644 --- a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp +++ b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp @@ -47,23 +47,31 @@ template struct PressureBoundaryRestrictedVelocityTo {}; } -template -struct component> { +template +struct component> { private: saw::data pressure_setting_; saw::data rho_setting_; public: component(const saw::data& pressure_setting__): pressure_setting_{pressure_setting__}, - rho_setting_{pressure_setting__ * df_info::inv_cs2} + rho_setting_{pressure_setting__ * df_info::inv_cs2} {} template - void apply(saw::data* field, saw::data> index, uint64_t time_step){ - using dfi = df_info; + void apply( + saw::data* field, + saw::data> index, + const saw::data>& meta, + uint64_t time_step + ) const { - bool is_even = ((time_step % 2) == 0); - auto& cell = field[index]; + using dfi = df_info; + + auto flat_ind = flatten_index::apply(index,meta); + + bool is_even = ((time_step % 2u) == 0u); + auto& cell = field[flat_ind.get()]; auto& info = cell.template get<"info">(); if(info({0u}).get() == 0u){ @@ -74,10 +82,11 @@ public: /** * Sum all known DFs */ - saw::data sum_df{0}; - for(saw::data k{0u}; k < saw::data{Descriptor::Q}; ++k){ + saw::data sum_df{0u}; + for(saw::data k{0u}; k < saw::data{Desc::Q}; ++k){ auto c_k = dfi::directions[k.get()]; - auto& cell_n = field({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}}); + auto flat_ind_n = flatten_index::apply({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}},meta); + auto& cell_n = field[flat_ind_n.get()]; auto& info_n = cell_n.template get<"info">(); auto info_n_val = info_n({0u}); auto k_opp = dfi::opposite_index[k.get()]; @@ -92,12 +101,12 @@ public: constexpr int known_dir = East ? 1 : -1; auto sum_unknown_dfs = (rho_setting_ - sum_df) * saw::data{known_dir}; - for(saw::data k{0u}; k < saw::data{Descriptor::Q}; ++k){ + for(saw::data k{0u}; k < saw::data{Desc::Q}; ++k){ auto c_k = dfi::directions[k.get()]; - auto& cell_n = field({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}}); + auto flat_ind_n = flatten_index::apply({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}},meta); + auto& cell_n = field[flat_ind_n.get()]; auto& info_n = cell_n.template get<"info">(); auto info_n_val = info_n({0u}); - // auto k_opp = dfi::opposite_index[k.get()]; if(info_n_val.get() > 0u){ sum_unknown_dfs += dfs_old({k}) * c_k[0u]; @@ -149,7 +158,6 @@ saw::error_or set_geometry( size_t j = idx[1]; size_t acc_id = j * meta.at({0u}).get() + i; - auto& c = cells[acc_id]; auto& info = c.template get<"info">()({0}); auto& dfs = c.template get<"dfs">(); @@ -165,7 +173,7 @@ saw::error_or set_geometry( // Left input info.set({3u}); }else if((i+2u) == meta.at({0u}).get() and (j >= 1 and (j+1) < meta.at({1u}).get() )){ - // Left output + // Right output info.set({4u}); }else { info.set({0u}); @@ -175,9 +183,7 @@ saw::error_or set_geometry( dfs_old(k) = eq.at(k); } }); - }); - - sycl_q.wait(); + }).wait(); return saw::make_void(); } @@ -193,7 +199,7 @@ void step( using namespace kel::lbm; using dfi = df_info; - constexpr saw::data frequency{1.0 / 0.59}; + constexpr saw::data frequency{1.0 / 0.5384}; bool is_even = ((time_step % 2u) == 0u); /** @@ -203,11 +209,15 @@ void step( component coll{0.5384}; component bb; */ - component> inlet{1.01 * dfi::cs2}; + component> inlet{1.1 * dfi::cs2}; component> outlet{1.0 * dfi::cs2}; + + // auto collision_ev = sycl_q.submit([&](acpp::sycl::handler& h){ - h.parallel_for(acpp::sycl::range<2>{meta.at({0}).get(), meta.at({1}).get()},[=](acpp::sycl::id<2> idx){ + + /// Collision + h.parallel_for(acpp::sycl::range<2>{meta.at({0}).get(), meta.at({1}).get()},[=](acpp::sycl::id idx){ size_t i = idx[0]; size_t j = idx[1]; size_t acc_id = j * meta.at({0u}).get() + i; @@ -216,18 +226,18 @@ void step( auto& info = cells[acc_id].template get<"info">(); switch (info({0u}).get()) { + // Bounce Back case 1u: { - // bb.apply(latt_acc, {i, j}, time_step); - auto& dfs_old = is_even ? c.template get<"dfs_old">() : c.template get<"dfs">(); auto df_cpy = dfs_old.copy(); - for(uint64_t i = 1u; i < Desc::Q; ++i){ - dfs_old({i}) = df_cpy({dfi::opposite_index.at(i)}); + for(uint64_t k = 1u; k < Desc::Q; ++k){ + dfs_old({k}) = df_cpy({dfi::opposite_index.at(k)}); } break; } + // Collision case 2u: { // coll.apply(latt_acc, {i, j}, time_step); auto& dfs_old = is_even ? c.template get<"dfs_old">() : c.template get<"dfs">(); @@ -236,23 +246,56 @@ void step( saw::data& rho = macro_c.template get<"pressure">(); saw::data>& vel = macro_c.template get<"velocity">(); + compute_rho_u(dfs_old,rho,vel); auto eq = equilibrium(rho,vel); - for(uint64_t i = 0u; i < Desc::Q; ++i){ - dfs_old({i}) = dfs_old({i}) + frequency * (eq.at(i) - dfs_old({i})); + for(uint64_t k = 0u; k < Desc::Q; ++k){ + dfs_old({k}) = dfs_old({k}) + frequency * (eq.at({k}) - dfs_old({k})); } break; } case 3u: { - inlet.apply( - }; + inlet.apply(cells, {{i,j}}, meta, time_step); + break; + } + case 4u: { + outlet.apply(cells, {{i,j}}, meta, time_step); + break; + } default: // Do nothing break; } }); - }); + }).wait(); + + //auto stream_ev = + sycl_q.submit([&](acpp::sycl::handler& h){ + /// Stream + h.parallel_for(acpp::sycl::range<2>{meta.at({0}).get(), meta.at({1}).get()},[=](acpp::sycl::id idx){ + size_t i = idx[0]; + size_t j = idx[1]; + size_t acc_id = j * meta.at({0u}).get() + i; + + auto& c = cells[acc_id]; + auto& info = c.template get<"info">(); + auto& dfs_new = is_even ? c.template get<"dfs">() : c.template get<"dfs_old">(); + + if (info({0u}).get() > 0u) { + for (uint64_t k = 0u; k < Desc::Q; ++k) { + auto dir = dfi::directions[dfi::opposite_index[k]]; + size_t acc_old_id = (j+dir[1]) * meta.at({0u}).get() + (i+dir[0]); + + auto& dfs_old = is_even ? cells[acc_old_id].template get<"dfs_old">() : cells[acc_old_id].template get<"dfs">(); + auto& info_old = cells[acc_old_id].template get<"info">(); + + dfs_new({k}) = dfs_old({k}); + } + } + + }); + }).wait(); } } } @@ -300,14 +343,14 @@ saw::error_or kel_main(int argc, char** argv){ } } - std::string vtk_f_name{"tmp/poiseulle_2d_gpu_"}; - vtk_f_name += std::to_string(0u) + ".vtk"; - // write_vtk_file(vtk_f_name,host_cells); for(uint64_t i = 0u; i < 1024u*1204u; ++i){ + sycl_q.wait(); if(i%128u == 0u){ - sycl_q.memcpy(&host_cells[0u], macro_cells, x_d * y_d * sizeof(saw::data) ); - sycl_q.wait(); + std::string vtk_f_name{"tmp/poiseulle_2d_gpu_"}; + vtk_f_name += std::to_string(i) + ".vtk"; + // write_vtk_file(vtk_f_name,host_cells); + sycl_q.memcpy(&host_cells[0u], macro_cells, x_d * y_d * sizeof(saw::data) ).wait(); lbm::write_vtk_file(vtk_f_name, &host_cells[0], meta); } lbm::step(cells,macro_cells,meta,i,sycl_q); -- cgit v1.2.3