From 869974982752cd37e808283af629ea0bbb680393 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Fri, 11 Jul 2025 14:14:36 +0200 Subject: Writing was flipped on non square, cubic domains --- examples/cavity_2d_gpu.cpp | 314 +++++++++++++++++++++++++++++++++++++++++++++ examples/poiseulle_2d.cpp | 222 +++++++++++++++++++++++++++++--- 2 files changed, 520 insertions(+), 16 deletions(-) create mode 100644 examples/cavity_2d_gpu.cpp (limited to 'examples') diff --git a/examples/cavity_2d_gpu.cpp b/examples/cavity_2d_gpu.cpp new file mode 100644 index 0000000..e1e6333 --- /dev/null +++ b/examples/cavity_2d_gpu.cpp @@ -0,0 +1,314 @@ +#include "../c++/descriptor.hpp" +#include "../c++/macroscopic.hpp" +#include "../c++/lbm.hpp" +#include "../c++/component.hpp" +#include "../c++/collision.hpp" +#include "../c++/boundary.hpp" + +#include +// #include +#include +#include + +namespace kel { +namespace lbm { +namespace sch { +using namespace saw::schema; + +/** + * Basic distribution function + * Base type + * D + * Q + * Scalar factor + * D factor + * Q factor + */ +using T = Float32; +using D2Q5 = Descriptor<2u,5u>; +using D2Q9 = Descriptor<2u,9u>; + +template +using DfCell = Cell; + +template +using CellInfo = Cell; + +/** + * Basic type for simulation + */ +template +using CellStruct = Struct< + Member, "dfs">, + Member, "dfs_old">, + Member, "info"> +>; + +template +using MacroStruct = Struct< + Member, "velocity">, + Member +>; + +using CavityFieldD2Q9 = CellField>; +} + + +/* +template +class df_cell_view; +*/ +/** + * Minor helper for the AA-Pull Pattern, so I can use only one lattice + * + * Am I sure I want to use AA this way? + * Esoteric Twist technically reduces the needed memory access footprint + */ +/* +template +class df_cell_view, Encode> { +public: + using Schema = sch::Cell; +private: + std::array::type>*, QN> view_; +public: + df_cell_view(const std::array::type>*, QN>& view): + view_{view} + {} +}; +*/ +namespace cmpt { +struct MovingWall {}; +} + +/** + * Full-Way moving wall Bounce back, something is not right here. + * Technically it should reflect properly. + */ +template +class component { +public: + std::array::type, Desc::D> lid_vel; + +public: + void apply( + saw::data>& dfs + ){ + using dfi = df_info; + + // Technically use .copy() + /* + auto dfs_cpy = dfs; + + for(uint64_t i = 0u; i < Desc::Q; ++i){ + dfs({dfi::opposite_index.at(i)}) = dfs_cpy({i}) - 2.0 * dfi::weights[i] * 1.0 * ( lid_vel[0] * dfi::directions[i][0] + lid_vel[1] * dfi::directions[i][1]) * dfi::inv_cs2; + } + */ + } +}; +} +} + +constexpr size_t dim_size = 2; +constexpr size_t dim_x = 128; +constexpr size_t dim_y = 128; + +void set_geometry(saw::data& latt){ + using namespace kel::lbm; + /** + * Set ghost + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& info = cell.template get<"info">(); + + info({0u}).set(0u); + + }, {{0u,0u}}, meta); + + /** + * Set wall + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& info = cell.template get<"info">(); + + info({0u}).set(2u); + + }, {{0u,0u}}, meta, {{1u,1u}}); + + /** + * Set fluid + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& info = cell.template get<"info">(); + + info({0u}).set(1u); + + }, {{0u,0u}}, meta, {{2u,2u}}); + + /** + * Set top lid + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& info = cell.template get<"info">(); + + info({0u}).set(3u); + + }, {{0u,1u}}, {{meta.at({0u}), 2u}}, {{2u,0u}}); +} + +void set_initial_conditions(saw::data& latt){ + using namespace kel::lbm; + + saw::data rho{1.0}; + saw::data> vel{{0.0,0.0}}; + auto eq = equilibrium(rho, vel); + + auto meta = latt.meta(); + + /** + * Set distribution + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& dfs = cell.template get<"dfs">(); + auto& dfs_old = cell.template get<"dfs_old">(); + + for(saw::data k = 0; k < saw::data{sch::D2Q9::Q}; ++k){ + dfs(k) = eq.at(k); + dfs_old(k) = eq.at(k); + } + + }, {{0u,0u}}, meta); +} + +void lbm_step( + saw::data& latt, + uint64_t time_step +){ + using namespace kel::lbm; + using dfi = df_info; + + /** + * 1. Relaxation parameter \tau + */ + component coll{0.5384}; + component bb; + + component bb_lid; + bb_lid.lid_vel = {0.1,0.0}; + + auto meta = latt.meta(); + + /** + * Collision + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& info = cell.template get<"info">(); + + auto& dfs = cell.template get<"dfs">(); + auto& dfs_old = cell.template get<"dfs_old">(); + + switch(info({0u}).get()){ + case 1u: { + coll.apply(latt, index, time_step); + break; + } + case 2u: { + bb.apply(latt, index, time_step); + break; + } + default: + break; + } + + }, {{0u,0u}}, meta); + + // Stream + for(uint64_t i = 1u; (i+1u) < latt.template get_dim_size<0>().get(); ++i){ + for(uint64_t j = 1u; (j+1u) < latt.template get_dim_size<1>().get(); ++j){ + auto& cell = latt({{i,j}}); + 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 && info_new({0u}).get() != 3u){ + for(uint64_t k = 0u; k < sch::D2Q9::Q; ++k){ + auto dir = dfi::directions[dfi::opposite_index[k]]; + auto& cell_dir_old = latt({{i+dir[0],j+dir[1]}}); + + auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">(); + auto& info_old = cell_dir_old.template get<"info">(); + + if( info_old({0}).get() == 3u ){ + auto& df_old_loc = even_step ? latt({{i,j}}).template get<"dfs_old">() : latt({{i,j}}).template get<"dfs">(); + df_new({k}) = df_old_loc({dfi::opposite_index.at(k)}) - 2.0 * dfi::inv_cs2 * dfi::weights.at(k) * 1.0 * ( bb_lid.lid_vel[0] * dir[0] + bb_lid.lid_vel[1] * dir[1]); + // dfs({dfi::opposite_index.at(i)}) = dfs_cpy({i}) - 2.0 * dfi::weights[i] * 1.0 * ( lid_vel[0] * dfi::directions[i][0] + lid_vel[1] * dfi::directions[i][1]) * dfi::inv_cs2; + } else { + df_new({k}) = df_old({k}); + } + } + } + } + } +} + +int main(){ + using namespace kel::lbm; + + saw::remote sycl_rmt; + + saw::data> dim{{dim_x, dim_y}}; + + saw::data lattice{dim}; + + converter conv{ + {0.1}, + {0.1} + }; + + print_lbm_meta(conv, {1e-3}); + + //auto& df_field = lattices.at(0).template get<"dfs">(); + //for(uint64_t i = 0; i < df_field.get_dim_size<0u>(); ++i){ + // lattices.at(i) = {dim_x, dim_y}; + //} + + /** + * Set meta information describing what this cell is + */ + set_geometry(lattice); + + /** + * + */ + set_initial_conditions(lattice); + + /** + * Timeloop + */ + + uint64_t lattice_steps = 512000u; + bool even_step = true; + + uint64_t print_every = 256u; + uint64_t file_no = 0u; + + saw::data,sch::D2Q9::D>> macros{dim}; + + for(uint64_t i = 0u; i < 256u; ++i){ + { + std::string vtk_f_name{"tmp/poiseulle_2d_"}; + vtk_f_name += std::to_string(i) + ".vtk"; + write_vtk_file(vtk_f_name, macros); + } + + // lbm_step(lattice, i); + } + return 0; +} diff --git a/examples/poiseulle_2d.cpp b/examples/poiseulle_2d.cpp index dcdc7c3..54623f7 100644 --- a/examples/poiseulle_2d.cpp +++ b/examples/poiseulle_2d.cpp @@ -45,8 +45,92 @@ using MacroStruct = Struct< Member >; +template +using GeometryStruct = Struct< + Member +>; + using CavityFieldD2Q9 = CellField>; } + +namespace cmpt { +template +struct PressureBoundaryRestrictedVelocityTo {}; +} + +/** + * This is massively hacky and expects a lot of conditions + * Either this or mirrored along the horizontal line works + * + * 0 - 2 - 2 + * 0 - 3 - 1 + * 0 - 3 - 1 + * ......... + * 0 - 3 - 1 + * 0 - 2 - 2 + * + */ +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__} + {} + + 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& 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">(); + + 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}); + } + } + constexpr int known_dir = East ? 1 : 1; + auto sum_unknown_dfs = (pressure_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 / pressure_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})); + } + } +}; } } @@ -114,8 +198,96 @@ void set_geometry(saw::data& latt){ void set_initial_conditions(saw::data& latt){ using namespace kel::lbm; + saw::data rho{1.0}; + saw::data> vel{{0.01,0.0}}; + auto eq = equilibrium(rho, vel); + auto meta = latt.meta(); + /** + * Set distribution + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& dfs = cell.template get<"dfs">(); + auto& dfs_old = cell.template get<"dfs_old">(); + + for(saw::data k = 0; k < saw::data{sch::D2Q9::Q}; ++k){ + dfs(k) = eq.at(k); + dfs_old(k) = eq.at(k); + } + + }, {{0u,0u}}, meta); +} + +void lbm_step( + saw::data& latt, + uint64_t time_step +){ + using namespace kel::lbm; + using dfi = df_info; + + bool even_step = ((time_step % 2u) == 0u); + /** + * 1. Relaxation parameter \tau + */ + component coll{0.5384}; + component bb; + component> inlet{1.1}; + component> outlet{0.9}; + + auto meta = latt.meta(); + + /** + * Collision + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& info = cell.template get<"info">(); + + auto& dfs = cell.template get<"dfs">(); + auto& dfs_old = cell.template get<"dfs_old">(); + + switch(info({0u}).get()){ + case 1u: { + coll.apply(latt, index, time_step); + break; + } + case 2u: { + bb.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 + iterate_over([&](const saw::data>& 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::D2Q9::Q; ++k){ + auto dir = dfi::directions[dfi::opposite_index[k]]; + auto& cell_dir_old = latt({{index.at({0u})+dir[0],index.at({1u})+dir[1]}}); + + 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); } int main(int argc, char** argv){ @@ -148,39 +320,57 @@ int main(int argc, char** argv){ print_lbm_meta>(conv, {conf.template get<"kinematic_viscosity">()}); - saw::data> dim{{32u, 8u}}; + saw::data> dim{{1024u, 128u}}; saw::data lattice{dim}; - - set_geometry(lattice); + auto meta = lattice.meta(); /** - * Print basic setup info + * Setup geometry */ - iterate_over([&](const saw::data>& index){ - auto& cell = lattice(index); - auto& info = cell.template get<"info">(); - auto info_v = info({0u}).get(); - std::cout<<(static_cast(info_v)); + set_geometry(lattice); - if( (index.at({1u}).get()+1u) < dim.at({1u}).get()){ - std::cout<<" "; - }else{ - std::cout<<"\n"; - } + { + + saw::data, sch::D2Q9::D>> geo{dim}; + + iterate_over([&](const saw::data>& index){ + auto& cell = lattice(index); + auto& info = cell.template get<"info">(); - }, {{0u,0u}}, dim); + geo(index).template get<"info">().set(info({0u}).get()); + }, {{0u,0u}}, dim); + + std::string vtk_f_name{"tmp/geometry.vtk"}; + write_vtk_file(vtk_f_name, geo); + } + + /** + * Setup DFs + */ + set_initial_conditions(lattice); saw::data,sch::D2Q9::D>> macros{dim}; for(uint64_t i = 0u; i < 256u; ++i){ + bool even_step = ((i % 2u) == 0u); { + // Stream + iterate_over([&](const saw::data>& index){ + auto& cell = lattice(index); + auto& dfs = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">(); + + auto& rho = macros.at(index).template get<"pressure">(); + auto& vel = macros.at(index).template get<"velocity">(); + compute_rho_u(dfs,rho,vel); + + }, {{0u,0u}}, meta); std::string vtk_f_name{"tmp/poiseulle_2d_"}; vtk_f_name += std::to_string(i) + ".vtk"; write_vtk_file(vtk_f_name, macros); } - // lbm_step(lattice, i); + lbm_step(lattice, i); } return 0; -- cgit v1.2.3