From 5a166a02efe011931faa7a9ff9fe4043e3a8017a Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Mon, 23 Jun 2025 16:12:27 +0200 Subject: Upgrading the cavity example --- examples/particle_ibm.cpp | 209 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 209 insertions(+) (limited to 'examples/particle_ibm.cpp') diff --git a/examples/particle_ibm.cpp b/examples/particle_ibm.cpp index a2b510e..0ce034e 100644 --- a/examples/particle_ibm.cpp +++ b/examples/particle_ibm.cpp @@ -1,7 +1,10 @@ #include "../c++/descriptor.hpp" +#include "../c++/equilibrium.hpp" #include +#include + namespace kel { namespace lbm { namespace sch { @@ -38,10 +41,196 @@ using CellStruct = Struct< Member >; +template +using MacroStruct = Struct< + Member, "velocity">, + Member +>; using CavityFieldD2Q9 = CellField>; } +namespace cmpt { +struct BounceBack{}; +struct MovingWall {}; +struct BGK {}; +struct ConstRhoBGK {}; } + +/** + * A reason for why a component based design is really good can be seen in my LR solver example + * + * Add Expression Templates and you're golden. + */ +template +class component; + +/** + * Full-Way BounceBack. I suspect that the moving wall requires half-way bounce back. + */ +template +class component { +public: + + void apply(saw::data>& dfs){ + using dfi = df_info; + + // Technically use .copy() + auto df_cpy = dfs; + + for(uint64_t i = 1u; i < Desc::Q; ++i){ + dfs({i}) = df_cpy({dfi::opposite_index.at(i)}); + } + } +}; + +/** + * 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; + } + */ + } +}; + +template +class component { +public: + typename saw::native_data_type::type relaxation_; +public: + void apply(saw::data>& dfs){ + typename saw::native_data_type::type rho; + std::array::type, Desc::D> vel; + compute_rho_u(dfs,rho,vel); + auto eq = equilibrium(rho,vel); + + for(uint64_t i = 0u; i < Desc::Q; ++i){ + dfs({i}).set(dfs({i}).get() + (1.0 / relaxation_) * (eq[i] - dfs({i}).get())); + } + } +}; + +template +class component { +public: + typename saw::native_data_type::type relaxation_; + typename saw::native_data_type::type rho_; +public: + void apply(saw::data>& dfs){ + std::array::type, Desc::D> vel; + compute_const_rho_u(dfs,rho_,vel); + auto eq = equilibrium(rho_,vel); + + for(uint64_t i = 0u; i < Desc::Q; ++i){ + dfs({i}).set(dfs({i}).get() + (1.0 / relaxation_) * (eq[i] - dfs({i}).get())); + } + } +}; + +} +} + +void set_geometry(saw::data& latt){ + using namespace kel::lbm; + /* + apply_for_cells([](auto& cell, std::size_t i, std::size_t j){ + uint8_t val = 0; + if(j == 1){ + val = 2u; + } + if(i == 1 || (i+2) == dim_x || (j+2) == dim_y){ + val = 3u; + } + if(i > 1 && (i+2) < dim_x && j > 1 && (j+2) < dim_y){ + val = 1u; + } + if(i == 0 || j == 0 || (i+1) == dim_x || (j+1) == dim_y){ + val = 0u; + } + cell.template get<"info">()(0u).set(val); + }, latt); + */ +} + +void lbm_step( + saw::data& latt, + bool even_step +){ + using namespace kel::lbm; + using dfi = df_info; + + component coll; + coll.relaxation_ = 0.5384; + component rho_coll; + rho_coll.relaxation_ = 0.5384; + rho_coll.rho_ = 1.0; + + component bb; + component bb_lid; + bb_lid.lid_vel = {0.1,0.0}; + + // Collide + apply_for_cells([&](auto& cell, std::size_t i, std::size_t j){ + auto& df = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + auto& info = cell.template get<"info">(); + + auto info_val = info({0u}).get(); + switch(info_val){ + case 1u: + coll.apply(df); + break; + case 2u: + // bb.apply(df); + bb_lid.apply(df); + break; + case 3u: + bb.apply(df); + break; + } + }, latt); + + // Stream + for(uint64_t i = 1; (i+1) < latt.template get_dim_size<0>().get(); ++i){ + for(uint64_t j = 1; (j+1) < 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() != 2u){ + 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() == 2u ){ + 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(int argc, char** argv){ @@ -78,9 +267,29 @@ int main(int argc, char** argv){ saw::data> dim{{128, 128}}; saw::data lattice{dim}; + lbm::particle_system> system; + /** + * Set meta information describing what this cell is + */ + set_geometry(lattice); + uint64_t lattice_steps{128u}; + uint64_t print_every{64u}; + bool even_step{true}; + + uint64_t file_num{0u}; + + saw::data,sch::D2Q9::D>> macros{dim}; + + for(uint64_t step{0u}; step < lattice_steps; ++step){ + + lbm_step(lattice, even_step); + + even_step = not even_step; + + } return 0; } -- cgit v1.2.3