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 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 314 insertions(+) create mode 100644 examples/cavity_2d_gpu.cpp (limited to 'examples/cavity_2d_gpu.cpp') 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; +} -- cgit v1.2.3