diff options
Diffstat (limited to 'examples/poiseulle_3d/poiseulle_3d.cpp')
| -rw-r--r-- | examples/poiseulle_3d/poiseulle_3d.cpp | 164 |
1 files changed, 152 insertions, 12 deletions
diff --git a/examples/poiseulle_3d/poiseulle_3d.cpp b/examples/poiseulle_3d/poiseulle_3d.cpp index 68e1bb7..9220a1f 100644 --- a/examples/poiseulle_3d/poiseulle_3d.cpp +++ b/examples/poiseulle_3d/poiseulle_3d.cpp @@ -2,6 +2,56 @@ #include <forstio/remote/sycl/sycl.hpp> +namespace saw { +template<typename Desc, typename CellT, typename Encode> +class data<kel::lbm::sch::CellField<Desc, CellT>, encode::Sycl<Encode>> final { +public: + using Schema = kel::lbm::sch::CellField<Desc,CellT>; + using MetaSchema = typename meta_schema<Schema>::MetaSchema; +private: + data<schema::Array<CellT,Desc::D>, encode::Sycl<Encode>> inner_; +public: + data() = delete; + data(const data<MetaSchema,Encode>& inner_meta__, acpp::sycl::queue& q): + inner_{inner_meta__,q} + {} + + const data<MetaSchema, Encode> meta() const { + return inner_.dims(); + } + + template<uint64_t i> + data<schema::UInt64,Encode> get_dim_size() const { + static_assert(i < Desc::D, "Not enough dimensions"); + return inner_.template get_dim_size<i>(); + } + + const data<CellT>& operator()(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index)const{ + return inner_.at(index); + } + + data<CellT>& operator()(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index){ + return inner_.at(index); + } + + const data<CellT>& at(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index)const{ + return inner_.at(index); + } + + data<CellT>& at(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index){ + return inner_.at(index); + } + + data<schema::UInt64,Encode> internal_size() const { + return inner_.internal_size(); + } + + data<CellT,Encode>* internal_data() { + return inner_.internal_data(); + } +}; +} + namespace kel { namespace lbm { namespace sch { @@ -21,8 +71,8 @@ using CellStruct = Struct< >; } -template<typename StructFieldSchema> -saw::error_or<void> set_geometry(saw::data<StructFieldSchema>& latt){ +template<typename StructFieldSchema, typename Encode> +saw::error_or<void> set_geometry(saw::data<StructFieldSchema, Encode>& latt){ using namespace kel::lbm; auto meta = latt.meta(); @@ -109,11 +159,96 @@ saw::error_or<void> set_geometry(saw::data<StructFieldSchema>& latt){ return saw::make_void(); } -saw::error_or<void> set_initial_conditions(saw::data<StructFieldSchema>& latt){ +template<typename StructFieldSchema, typename Encode> +saw::error_or<void> set_initial_conditions(saw::data<StructFieldSchema, Encode>& latt){ + + saw::data<sch::T> rho{1.0}; + saw::data<sch::FixedArray<sch::T,sch::Desc::D>> vel{}; + vel.at({0u}) = {0.1f}; + auto eq = equilibrium<sch::T,sch::Desc>(rho, vel); + + auto meta = latt.meta(); + + /** + * Set distribution + */ + iterator<sch::Desc::D>::apply([&](const saw::data<sch::FixedArray<sch::UInt64,sch::Desc::D>>& index){ + auto& cell = latt.at(index); + auto& dfs = cell.template get<"dfs">(); + auto& dfs_old = cell.template get<"dfs_old">(); + + for(saw::data<sch::UInt64> k = 0; k < saw::data<sch::UInt64>{sch::Desc::Q}; ++k){ + dfs(k) = eq.at(k); + dfs_old(k) = eq.at(k); + } + }, {}, meta); return saw::make_void(); } -saw::error_or<void> lbm_step(saw::data<StructFieldSchema>& latt){ +template<typename StructFieldSchema, typename Encode> +saw::error_or<void> lbm_step( + saw::data<StructFieldSchema, Encode>& latt, + saw::data<sch::UInt64> time_step, + acpp::sycl::queue& sycl_q +){ + using dfi = df_info<sch::T,sch::Desc>; + auto meta = latt.meta(); + bool even_step = ((time_step.get() % 2u) == 0u); + + component<sch::T, sch::Desc, cmpt::BGK, saw::encode::Sycl<saw::encode::Native>> coll{0.5384}; + component<sch::T, sch::Desc, cmpt::BounceBack, saw::encode::Sycl<saw::encode::Native>> bb; + + /** + * Collision + */ + iterator<sch::Desc::D>::apply([&](const saw::data<sch::FixedArray<sch::UInt64,sch::Desc::D>>& index){ + auto& cell = latt.at(index); + auto& info = cell.template get<"info">(); + + switch(info({0u}).get()){ + case 1u: { + bb.apply(latt, index, time_step); + break; + } + case 2u: { + coll.apply(latt, index, time_step); + break; + } + case 3u: { + //inlet.apply(latt, index, time_step); + break; + } + case 4u: { + //outlet.apply(latt, index, time_step); + break; + } + default: + break; + } + + }, {{0u,0u}}, meta); + + // Stream + iterator<sch::Desc::D>::apply([&](const saw::data<sch::FixedArray<sch::UInt64,sch::Desc::D>>& index){ + auto& cell = latt(index); + 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& cell_dir_old = latt.at(index_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}); + } + } + }, {{0u,0u}}, meta); return saw::make_void(); } @@ -126,31 +261,36 @@ saw::error_or<void> real_main(int argc, char** argv){ }; print_lbm_meta<sch::T,sch::Desc>(conv, {1e-3f}); + acpp::sycl::queue sycl_q; /** * Set the meta and initialize the lattice */ saw::data<sch::FixedArray<sch::UInt64,sch::Desc::D>> meta{{2048u,128u,128u}}; - saw::data<sch::CellField<sch::Desc,sch::CellStruct>, encode::Sycl<encode::Native>> lattice{meta}; + saw::data<sch::CellField<sch::Desc,sch::CellStruct>, saw::encode::Sycl<saw::encode::Native>> lattice{meta,sycl_q}; /** * Set the geometry */ - auto eov = set_geometry(lattice); - if(eov.is_error()){ - return eov; + { + auto eov = set_geometry(lattice); + if(eov.is_error()){ + return eov; + } } /** * Set initial conditions in the simulations */ + { auto eov = set_initial_conditions(lattice); - if(eov.is_error()){ - return eov; + if(eov.is_error()){ + return eov; + } } - for(uint64_t i = 0u; i < 1024u*32; ++i){ - auto eov = lbm_step(lattice); + for(saw::data<sch::UInt64> ts{0u}; ts < saw::data<sch::UInt64>{1024u*32u}; ++ts){ + auto eov = lbm_step(lattice,ts,sycl_q); if(eov.is_error()){ return eov; } |