#include #include namespace saw { /* template class data, encode::Sycl> final { public: using Schema = kel::lbm::sch::CellField; using MetaSchema = typename meta_schema::MetaSchema; private: data, encode::Sycl> inner_; public: data() = delete; data(const data& inner_meta__, acpp::sycl::queue& q): inner_{inner_meta__,q} {} const data meta() const { return inner_.dims(); } template data get_dim_size() const { static_assert(i < Desc::D, "Not enough dimensions"); return inner_.template get_dim_size(); } const data& operator()(const data, Encode>& index)const{ return inner_.at(index); } data& operator()(const data, Encode>& index){ return inner_.at(index); } const data& at(const data, Encode>& index)const{ return inner_.at(index); } data& at(const data, Encode>& index){ return inner_.at(index); } ref, encode::Sycl >> inner_data() { return {inner_}; } data internal_size() const { return inner_.internal_size(); } data* internal_data() { return inner_.internal_data(); } }; */ /* template class data>, encode::Sycl> final { public: using Schema = schema::Ref>; private: data>, encode::Sycl>> inner_ref_; public: data() = delete; data(ref,encode::Sycl>> inner_ref__): inner_ref_{inner_ref__().inner_data()} {} const data meta() const { return inner_ref_().dims(); } const data& operator()(const data, Encode>& index)const{ return inner_ref_().at(index); } data& operator()(const data, Encode>& index){ return inner_ref_().at(index); } const data& at(const data, Encode>& index)const{ return inner_ref_().at(index); } data& at(const data, Encode>& index){ return inner_ref_().at(index); } data internal_size() const { return inner_ref_().internal_size(); } data* internal_data() { return inner_ref_().internal_data(); } }; template class data>, encode::Sycl> final { public: using Schema = schema::Ref>; private: data>, encode::Sycl>> inner_ref_; public: data() = delete; data(ref,encode::Sycl>> inner_ref__): inner_ref_{inner_ref__().inner_data()} {} const data meta() const { return inner_ref_().dims(); } const data& operator()(const data, Encode>& index)const{ return inner_ref_().at(index); } data& operator()(const data, Encode>& index){ return inner_ref_().at(index); } const data& at(const data, Encode>& index)const{ return inner_ref_().at(index); } data& at(const data, Encode>& index){ return inner_ref_().at(index); } data internal_size() const { return inner_ref_().internal_size(); } data* internal_data() { return inner_ref_().internal_data(); } }; */ } namespace kel { namespace lbm { namespace sch { using namespace saw::schema; using T = Float32; using Desc = Descriptor<3u,27u>; using DfCell = Cell; using CellInfo = Cell; using CellStruct = Struct< Member, Member, //Member, "force">, Member >; } using SyclKelAllocator = acpp::sycl::usm_allocator, acpp::sycl::usm::alloc::shared>; uint64_t to_flat_index(const saw::data>& index, const saw::data>& meta){ auto flat = index.at({2u}) + meta.at({2u})*(index.at({1u}) + index.at({0u})*meta.at({1u})); return flat.get(); } template saw::error_or set_geometry(std::vector,SyclKelAllocator>& latt, const saw::data>& meta){ using namespace kel::lbm; /** * Set ghost */ iterator::apply([&](const saw::data>& index){ auto findex = to_flat_index(index,meta); auto& cell = latt.at(findex); auto& info = cell.template get<"info">(); info({0u}).set(0u); }, {{0u,0u}}, meta); saw::data> distance; /** * Set wall */ iterator<1u>::apply([&](const saw::data>& index){ distance.at(index.at({0u})) = {1u}; },{{0u}},{{sch::Desc::D}}); /* */ iterator::apply([&](const saw::data>& index){ auto findex = to_flat_index(index,meta); auto& cell = latt.at(findex); auto& info = cell.template get<"info">(); info({0u}).set(2u); }, {}, meta, distance); /** * Set fluid */ iterator<1u>::apply([&](const saw::data>& index){ distance.at(index.at({0u})) = {2u}; },{{0u}},{{sch::Desc::D}}); /* */ iterator::apply([&](const saw::data>& index){ auto findex = to_flat_index(index,meta); auto& cell = latt.at(findex); auto& info = cell.template get<"info">(); info({0u}).set(1u); }, {}, meta, distance); /** * Set inflow */ saw::data> shifted_point; iterator<1u>::apply([&](const saw::data>& index){ distance.at(index.at({0u})) = (index.at({0u}).get() == 0u)?1u:2u; shifted_point.at(index.at({0u})) = (index.at({0u}).get() == 0u)?3u:meta.at({index.at({0u})}); },{{0u}},{{sch::Desc::D}}); /* */ iterator::apply([&](const saw::data>& index){ auto findex = to_flat_index(index,meta); auto& cell = latt.at(findex); auto& info = cell.template get<"info">(); info({0u}).set(3u); }, {}, shifted_point, distance); /** * Set outflow */ iterator<1u>::apply([&](const saw::data>& index){ distance.at(index.at({0u})) = (index.at({0u}).get() == 0u)?1u:2u; shifted_point.at(index.at({0u})) = (index.at({0u}).get() == 0u)?(meta.at({0u})-3u):0u; },{{0u}},{{sch::Desc::D}}); /* */ iterator::apply([&](const saw::data>& index){ auto findex = to_flat_index(index,meta); auto& cell = latt.at(findex); auto& info = cell.template get<"info">(); info({0u}).set(4u); }, shifted_point, meta, distance); return saw::make_void(); } template saw::error_or set_initial_conditions(std::vector,SyclKelAllocator>& latt, const saw::data>& meta){ saw::data rho{1.0}; saw::data> vel{}; vel.at({0u}) = {0.1f}; auto eq = equilibrium(rho, vel); /** * Set distribution */ iterator::apply([&](const saw::data>& index){ auto findex = to_flat_index(index,meta); auto& cell = latt.at(findex); auto& dfs = cell.template get<"dfs">(); auto& dfs_old = cell.template get<"dfs_old">(); for(saw::data k = 0; k < saw::data{sch::Desc::Q}; ++k){ dfs(k) = eq.at(k); dfs_old(k) = eq.at(k); } }, {}, meta); return saw::make_void(); } template saw::error_or lbm_step( std::vector,SyclKelAllocator>& latt, const saw::data>& meta, saw::data time_step, acpp::sycl::queue& sycl_q ){ using dfi = df_info; bool even_step = ((time_step.get() % 2u) == 0u); component coll{0.5384}; component bb; /** * Collision */ iterator::apply([&](const saw::data>& index){ auto findex = to_flat_index(index,meta); auto& cell = latt.at(findex); auto& info = cell.template get<"info">(); switch(info({0u}).get()){ case 1u: { bb.apply(&latt[0u], meta, index, time_step); break; } case 2u: { coll.apply(&latt[0u], meta, index, time_step); break; } case 3u: { //inlet.apply(latt, index, time_step); break; } case 4u: { //outlet.apply(latt, index, time_step); break; } default: break; } }, {}, meta); // Stream iterator::apply([&](const saw::data>& index){ auto findex = to_flat_index(index,meta); auto& cell = latt.at(findex); auto& df_new = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">(); auto& info_new = cell.template get<"info">(); if(info_new({0u}).get() > 0u){ for(uint64_t k = 0u; k < sch::Desc::Q; ++k){ auto dir = dfi::directions[dfi::opposite_index[k]]; auto index_dir = index; for(uint64_t i = 0; i < sch::Desc::D; ++i){ index_dir.at({i}) = index_dir.at({i}) + dir[i]; } auto findex_dir = to_flat_index(index_dir,meta); auto& cell_dir_old = latt.at(findex_dir); auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">(); df_new({k}) = df_old({k}); } } }, {}, meta); return saw::make_void(); } saw::error_or real_main(int argc, char** argv){ using dfi = df_info; converter conv { {0.1f}, {0.1f} }; print_lbm_meta(conv, {1e-3f}); acpp::sycl::queue sycl_q; /** * Set the meta and initialize the lattice */ saw::data> meta{{2048u,128u,128u}}; SyclKelAllocator sycl_shared_alloc{sycl_q}; std::vector,SyclKelAllocator> lattice{(meta.at({0u})*meta.at({1u})*meta.at({2u})).get(), sycl_shared_alloc}; /** * Set the geometry */ { auto eov = set_geometry(lattice, meta); if(eov.is_error()){ return eov; } } /** * Set initial conditions in the simulations */ { auto eov = set_initial_conditions(lattice, meta); if(eov.is_error()){ return eov; } } for(saw::data ts{0u}; ts < saw::data{1024u*32u}; ++ts){ auto eov = lbm_step(lattice,meta,ts,sycl_q); if(eov.is_error()){ return eov; } } return saw::make_void(); } } // lbm } // kel /** * main, but I don't like the error handling */ int main(int argc, char** argv){ auto eov = kel::lbm::real_main(argc, argv); if(eov.is_error()){ auto& err = eov.get_error(); auto err_msg = err.get_message(); std::cerr<<"[Error]: "<