From 6f6ea556a54f46716877ae44adcd456aec56226e Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Fri, 19 Dec 2025 20:18:09 +0100 Subject: Adding random stuff --- examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp | 172 +++++++++++++++++++++---- 1 file changed, 148 insertions(+), 24 deletions(-) (limited to 'examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp') diff --git a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp index 6dda469..b31f8c5 100644 --- a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp +++ b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp @@ -6,9 +6,6 @@ template using SyclHostAlloc = acpp::sycl::usm_allocator, acpp::sycl::usm::alloc::host>; -template -using SyclDeviceAlloc = acpp::sycl::usm_allocator, acpp::sycl::usm::alloc::device>; - namespace kel { namespace lbm { namespace sch { @@ -38,10 +35,9 @@ using CellStruct = Struct< Member >; -template using MacroStruct = Struct< - Member, "velocity">, - Member + Member, "velocity">, + Member >; } @@ -51,6 +47,77 @@ template struct PressureBoundaryRestrictedVelocityTo {}; } +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} + {} + + template + void apply(saw::data* field, saw::data> index, uint64_t time_step){ + using dfi = df_info; + + bool is_even = ((time_step % 2) == 0); + auto& cell = field[index]; + + auto& info = cell.template get<"info">(); + if(info({0u}).get() == 0u){ + return; + } + auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + + /** + * Sum all known DFs + */ + saw::data sum_df{0}; + for(saw::data k{0u}; k < saw::data{Descriptor::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& 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_df += dfs_old({k_opp}); + } + } + /** + * Get the sum of the unknown dfs and precalculate the direction + */ + 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){ + auto c_k = dfi::directions[k.get()]; + auto& cell_n = field({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}}); + 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]; + } + } + + auto vel_x = sum_unknown_dfs / rho_setting_; + + if constexpr (East) { + dfs_old({2u}) = dfs_old({1u}) + saw::data{2.0 / 3.0} * rho_setting_ * vel_x; + dfs_old({6u}) = dfs_old({5u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({4u}) - dfs_old({3u})); + dfs_old({8u}) = dfs_old({7u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({3u}) - dfs_old({4u})); + }else if constexpr (not East){ + dfs_old({1u}) = dfs_old({2u}) - saw::data{2.0 / 3.0} * rho_setting_ * vel_x; + dfs_old({5u}) = dfs_old({6u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({3u}) - dfs_old({4u})); + dfs_old({7u}) = dfs_old({8u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({4u}) - dfs_old({3u})); + } + } +}; + /** * This is massively hacky and expects a lot of conditions * Either this or mirrored along the horizontal line works @@ -64,23 +131,25 @@ struct PressureBoundaryRestrictedVelocityTo {}; * */ +template saw::error_or set_geometry( saw::data* cells, - const saw::data>& meta, + const saw::data>& meta, acpp::sycl::queue& sycl_q ){ using namespace kel::lbm; saw::data rho{1.0}; - saw::data> vel{{0.0,0.0}}; - auto eq = equilibrium(rho, vel); + saw::data> vel{{0.0,0.0}}; + auto eq = equilibrium(rho, vel); 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){ 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">()({0}); auto& dfs = c.template get<"dfs">(); @@ -113,25 +182,77 @@ saw::error_or set_geometry( return saw::make_void(); } -void lbm_step( +template +void step( saw::data* cells, - const saw::data>& meta, + saw::data* macro_cells, + const saw::data>& meta, uint64_t time_step, acpp::sycl::queue& sycl_q ){ using namespace kel::lbm; - using dfi = df_info; + using dfi = df_info; + + constexpr saw::data frequency{1.0 / 0.59}; - bool even_step = ((time_step % 2u) == 0u); + bool is_even = ((time_step % 2u) == 0u); /** * 1. Relaxation parameter \tau */ /* component coll{0.5384}; component bb; - component> inlet{1.01 * dfi::cs2}; - component> outlet{1.0 * dfi::cs2}; */ + component> inlet{1.01 * dfi::cs2}; + component> outlet{1.0 * dfi::cs2}; + + 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){ + 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 = cells[acc_id].template get<"info">(); + + switch (info({0u}).get()) { + 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)}); + } + + break; + } + case 2u: { + // coll.apply(latt_acc, {i, j}, time_step); + auto& dfs_old = is_even ? c.template get<"dfs_old">() : c.template get<"dfs">(); + + auto& macro_c = macro_cells[acc_id]; + + 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})); + } + break; + } + case 3u: { + inlet.apply( + }; + default: + // Do nothing + break; + } + }); + }); } } } @@ -162,16 +283,18 @@ saw::error_or kel_main(int argc, char** argv){ uint64_t x_d = 256u; uint64_t y_d = 64u; + saw::data> meta{{x_d,y_d}}; + acpp::sycl::queue sycl_q; - SyclHostAlloc sycl_host_alloc{sycl_q}; + SyclHostAlloc sycl_host_alloc{sycl_q}; // SyclDeviceAlloc sycl_dev_alloc{sycl_q}; - std::vector, SyclHostAlloc> host_cells{x_d * y_d,sycl_host_alloc}; + std::vector, SyclHostAlloc> host_cells{x_d * y_d,sycl_host_alloc}; saw::data* cells = acpp::sycl::malloc_device>(x_d * y_d,sycl_q); - + saw::data* macro_cells = acpp::sycl::malloc_device>(x_d * y_d,sycl_q); { - auto eov = lbm::set_geometry(cells,{{x_d,y_d}},sycl_q); + auto eov = lbm::set_geometry(cells,meta,sycl_q); if(eov.is_error()){ return eov; } @@ -183,16 +306,17 @@ saw::error_or kel_main(int argc, char** argv){ for(uint64_t i = 0u; i < 1024u*1204u; ++i){ if(i%128u == 0u){ - + sycl_q.memcpy(&host_cells[0u], macro_cells, x_d * y_d * sizeof(saw::data) ); + sycl_q.wait(); + lbm::write_vtk_file(vtk_f_name, &host_cells[0], meta); } - lbm_step(cells,meta,i,sycl_q); + lbm::step(cells,macro_cells,meta,i,sycl_q); } - sycl_q.wait(); - sycl_q.memcpy(&host_cells[0u], cells, x_d * y_d * sizeof(saw::data) ); sycl_q.wait(); acpp::sycl::free(cells, sycl_q); + acpp::sycl::free(macro_cells, sycl_q); sycl_q.wait(); return saw::make_void(); -- cgit v1.2.3