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 --- c++/component.hpp | 9 +- examples/SConscript | 8 +- examples/cavity_2d.cpp | 43 ++------ examples/particle_ibm.cpp | 209 +++++++++++++++++++++++++++++++++++++++ tests/particle_flow_coupling.cpp | 20 ++++ 5 files changed, 252 insertions(+), 37 deletions(-) create mode 100644 tests/particle_flow_coupling.cpp diff --git a/c++/component.hpp b/c++/component.hpp index c67387b..87ac31c 100644 --- a/c++/component.hpp +++ b/c++/component.hpp @@ -6,6 +6,10 @@ namespace kel { namespace lbm { namespace cmpt { + +/** + * Maybe for the future + */ template class access { public: @@ -16,6 +20,9 @@ public: static constexpr bool skip_sync = SkipSync; }; +/** + * Maybe for the future + */ template class access_tuple { public: @@ -23,7 +30,7 @@ public: } /** - * Compponent class which forms the basis of the + * Component class which forms the basis of each processing step */ template class component {}; diff --git a/examples/SConscript b/examples/SConscript index 73adcf5..13ed580 100644 --- a/examples/SConscript +++ b/examples/SConscript @@ -28,10 +28,16 @@ examples_objects = []; examples_env.add_source_files(examples_objects, ['cavity_2d.cpp'], shared=False); examples_env.cavity_2d = examples_env.Program('#bin/cavity_2d', [env.library_static, examples_objects]); +# Channel Throat 2D +#examples_objects = []; +#examples_env.add_source_files(examples_objects, ['particle_ibm.cpp'], shared=False); +#examples_env.particle_ibm_2d = examples_env.Program('#bin/particle_ibm_2d', [env.library_static, examples_objects]); + # Set Alias env.examples = [ examples_env.meta_2d, - examples_env.cavity_2d + examples_env.cavity_2d, +# examples_env.particle_ibm_2d ]; env.Alias('examples', env.examples); diff --git a/examples/cavity_2d.cpp b/examples/cavity_2d.cpp index a10276a..654a9ac 100644 --- a/examples/cavity_2d.cpp +++ b/examples/cavity_2d.cpp @@ -1,6 +1,7 @@ #include "../c++/descriptor.hpp" #include "../c++/macroscopic.hpp" #include "../c++/lbm.hpp" +#include "../c++/component.hpp" /** */ @@ -79,14 +80,6 @@ std::array::type,Desc::Q> equilibrium( return eq; } -/** - * 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; - /* template class df_cell_view; @@ -122,7 +115,7 @@ struct ConstRhoBGK {}; * Full-Way BounceBack. I suspect that the moving wall requires half-way bounce back. */ template -class component { +class component { public: void apply(saw::data>& dfs){ @@ -142,7 +135,7 @@ public: * Technically it should reflect properly. */ template -class component { +class component { public: std::array::type, Desc::D> lid_vel; @@ -164,7 +157,7 @@ public: }; template -class component { +class component { public: typename saw::native_data_type::type relaxation_; public: @@ -179,23 +172,6 @@ public: } } }; - -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())); - } - } -}; } } @@ -284,14 +260,11 @@ void lbm_step( using namespace kel::lbm; using dfi = df_info; - component coll; + component coll; coll.relaxation_ = 0.5384; - component rho_coll; - rho_coll.relaxation_ = 0.5384; - rho_coll.rho_ = 1.0; - component bb; - component bb_lid; + component bb; + component bb_lid; bb_lid.lid_vel = {0.1,0.0}; // Collide @@ -505,7 +478,7 @@ int main(){ lbm_step(lattice, even_step); - even_step = !even_step; + even_step = not even_step; } /** 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; } diff --git a/tests/particle_flow_coupling.cpp b/tests/particle_flow_coupling.cpp new file mode 100644 index 0000000..216c229 --- /dev/null +++ b/tests/particle_flow_coupling.cpp @@ -0,0 +1,20 @@ +#include + +#include +#include "../c++/particle/geometry/circle.hpp" + + +namespace { +namespace sch { +using namespace kel::lbm::sch; + +using T = Float64; +} + +SAW_TEST("Particle Coupling"){ + using namespace kel; + lbm::particle_system> system; + + +} +} -- cgit v1.2.3