From a7420c5f5f56bb21de0241ed152ad9b55965d42d Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Fri, 30 Jan 2026 19:08:54 +0100 Subject: Simulation works partially. Reconfirm with writing data out --- examples/poiseulle_particles_2d_gpu/sim.cpp | 52 +++++++++++++++-- lib/core/c++/boundary.hpp | 90 +++++++++-------------------- lib/core/c++/chunk.hpp | 13 +++-- lib/core/c++/collision.hpp | 1 + lib/core/c++/lbm.hpp | 1 + lib/core/c++/stream.hpp | 25 ++++++++ 6 files changed, 109 insertions(+), 73 deletions(-) create mode 100644 lib/core/c++/stream.hpp diff --git a/examples/poiseulle_particles_2d_gpu/sim.cpp b/examples/poiseulle_particles_2d_gpu/sim.cpp index fc9fc61..e655315 100644 --- a/examples/poiseulle_particles_2d_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_gpu/sim.cpp @@ -32,7 +32,7 @@ saw::error_or setup_initial_conditions(saw::data> // Set everything as walls iterator::apply( [&](auto& index){ - info_f.at(index).set(2u); + info_f.at(index).set(1u); }, {}, info_f.get_dims(), @@ -41,7 +41,7 @@ saw::error_or setup_initial_conditions(saw::data> // Fluid iterator::apply( [&](auto& index){ - info_f.at(index).set(1u); + info_f.at(index).set(2u); }, {}, info_f.get_dims(), @@ -91,15 +91,59 @@ saw::error_or step(saw::data>,encode::Sy auto& q = dev.get_handle(); auto& info_f = fields.template get<"info">(); + // auto coll_ev = q.submit([&](acpp::sycl::handler& h){ + // Need nicer things to handle the flow. I see improvement here component> collision{0.6}; + component> bb; + component,encode::Sycl> flow_in{1.01}; + component,encode::Sycl> flow_out{1.0}; + h.parallel_for(acpp::sycl::range{dim_x,dim_y}, [=](acpp::sycl::id idx){ saw::data> index; for(uint64_t i = 0u; i < Desc::D; ++i){ index.at({{i}}).set(idx[i]); } - collision.apply(fields,index,t_i); + auto info = info_f.at(index); + + switch(info.get()){ + case 0u: + break; + case 1u: + bb.apply(fields,index,t_i); + break; + case 2u: + collision.apply(fields,index,t_i); + break; + case 3u: + flow_in.apply(fields,index,t_i); + collision.apply(fields,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + collision.apply(fields,index,t_i); + break; + default: + break; + } + }); + }).wait(); + + q.submit([&](acpp::sycl::handler& h){ + component> stream; + + h.parallel_for(acpp::sycl::range{dim_x,dim_y}, [=](acpp::sycl::id idx){ + saw::data> index; + for(uint64_t i = 0u; i < Desc::D; ++i){ + index.at({{i}}).set(idx[i]); + } + + auto info = info_f.at(index); + + if(info.get() > 0u){ + stream.apply(fields,index,t_i); + } }); }).wait(); @@ -155,7 +199,7 @@ saw::error_or lbm_main(int argc, char** argv){ } sycl_q.wait(); - for(saw::data i{0u}; i < saw::data{32ul}; ++i){ + for(saw::data i{0u}; i < saw::data{1024ul}; ++i){ auto eov = step(lsd_view,i,dev); if(eov.is_error()){ return eov; diff --git a/lib/core/c++/boundary.hpp b/lib/core/c++/boundary.hpp index 8ab2457..ec16edf 100644 --- a/lib/core/c++/boundary.hpp +++ b/lib/core/c++/boundary.hpp @@ -44,47 +44,22 @@ public: * Raw setup */ template - void apply(saw::data& field, const saw::data>& index, saw::data time_step){ + void apply(const saw::data& field, const saw::data>& index, saw::data time_step)const{ bool is_even = ((time_step.get() % 2u) == 0u); // This is a ref - auto& cell = field(index); - auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); - // auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + // auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">(); + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); using dfi = df_info; // Technically use .copy() - auto df_cpy = dfs_old.copy(); + auto& dfs_old = dfs_old_f.at(index); + auto df_cpy = dfs_old; - for(uint64_t i = 1u; i < Descriptor::Q; ++i){ - dfs_old({i}) = df_cpy({dfi::opposite_index.at(i)}); - } - } - - template - void apply( - saw::data* field, - const saw::data>& meta, - const saw::data>& index, - saw::data time_step - ){ - bool is_even = ((time_step.get() % 2u) == 0u); - - // This is a ref - auto& cell = field[0u]; - - auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); - // auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); - - using dfi = df_info; - - // Technically use .copy() - auto df_cpy = dfs_old.copy(); - - for(uint64_t i = 1u; i < Descriptor::Q; ++i){ - dfs_old({i}) = df_cpy({dfi::opposite_index.at(i)}); + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + dfs_old.at({i}) = df_cpy.at({dfi::opposite_index.at(i)}); } } }; @@ -104,59 +79,46 @@ public: template class component, Encode> final { 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} + component(const saw::data& density_setting__): + rho_setting_{density_setting__} {} template - void apply(saw::data& field, saw::data> index, saw::data time_step){ + void apply(const saw::data& field, saw::data> index, saw::data time_step) const { using dfi = df_info; + constexpr int known_dir = Dir ? 1 : -1; bool is_even = ((time_step.get() % 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">(); - // auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + auto& info_f = field.template get<"info">(); + + // auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">(); + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + auto& dfs_old = dfs_old_f.at(index); /** * 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}); + if(c_k[0u]*known_dir <= 0){ + sum_df += dfs_old.at({k}); } } /** * Get the sum of the unknown dfs and precalculate the direction */ - constexpr int known_dir = Dir ? 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]; + if(c_k[0u]*known_dir > 0){ + sum_unknown_dfs += dfs_old.at({k}) * c_k[0u]; } } @@ -165,13 +127,13 @@ public: static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly"); if constexpr (Dir) { - 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})); + dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data{2.0 / 3.0} * rho_setting_ * vel_x; + dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); + dfs_old.at({8u}) = dfs_old.at({7u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); }else if constexpr (not Dir){ - 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})); + dfs_old.at({1u}) = dfs_old.at({2u}) - saw::data{2.0 / 3.0} * rho_setting_ * vel_x; + dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); + dfs_old.at({7u}) = dfs_old.at({8u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); } } }; diff --git a/lib/core/c++/chunk.hpp b/lib/core/c++/chunk.hpp index 5d20faa..a1f2451 100644 --- a/lib/core/c++/chunk.hpp +++ b/lib/core/c++/chunk.hpp @@ -35,6 +35,7 @@ using SuperChunk = Array; } namespace saw { + template class data,Encode> final { public: @@ -57,19 +58,21 @@ public: return data::get_dims(); } - data& at(const data>& index){ + static constexpr auto to_ghost_index(const data>& index){ std::decay_t ind; for(uint64_t i = 0u; i < sizeof...(Sides); ++i){ ind.at({i}) = index.at({i}) + Ghost; } + return ind; + } + + data& at(const data>& index){ + std::decay_t ind = to_ghost_index(index); return values_.at(ind); } const data& at(const data>& index) const { - std::decay_t ind; - for(uint64_t i = 0u; i < sizeof...(Sides); ++i){ - ind.at({i}) = index.at({i}) + Ghost; - } + std::decay_t ind = to_ghost_index(index); return values_.at(ind); } diff --git a/lib/core/c++/collision.hpp b/lib/core/c++/collision.hpp index 85c3af9..ed26a08 100644 --- a/lib/core/c++/collision.hpp +++ b/lib/core/c++/collision.hpp @@ -3,6 +3,7 @@ #include "macroscopic.hpp" #include "component.hpp" #include "equilibrium.hpp" +#include "chunk.hpp" namespace kel { namespace lbm { diff --git a/lib/core/c++/lbm.hpp b/lib/core/c++/lbm.hpp index 00f153a..5334e4e 100644 --- a/lib/core/c++/lbm.hpp +++ b/lib/core/c++/lbm.hpp @@ -13,6 +13,7 @@ #include "equilibrium.hpp" #include "iterator.hpp" #include "macroscopic.hpp" +#include "stream.hpp" #include "write_vtk.hpp" #include "util.hpp" diff --git a/lib/core/c++/stream.hpp b/lib/core/c++/stream.hpp new file mode 100644 index 0000000..d217373 --- /dev/null +++ b/lib/core/c++/stream.hpp @@ -0,0 +1,25 @@ +#pragma once + +#include "component.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +struct Stream {}; +} + +template +class component final { +private: +public: + static constexpr saw::string_literal name = "streaming"; + static constexpr saw::string_literal after = "collide"; + static constexpr saw::string_literal before = ""; + + template + void apply(const saw::data& field, const saw::data>& index, saw::data time_step) const { + + } +}; +} +} -- cgit v1.2.3