From 04940fe64b714bcfad9026d54ee75a2fbdc7fe7b Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Wed, 17 Sep 2025 21:35:41 +0200 Subject: Streamlining and thoughts about changes in particle dynamics --- c++/descriptor.hpp | 1 + c++/macroscopic.hpp | 29 +++++++++++++++++++++ c++/util.hpp | 40 ++++++++++++++++++++++------- c++/write_vtk.hpp | 26 +++++++++++++++++++ examples/poiseulle_particles_channel_2d.cpp | 29 +++++++++++++++------ 5 files changed, 108 insertions(+), 17 deletions(-) diff --git a/c++/descriptor.hpp b/c++/descriptor.hpp index 0c6707a..d04d217 100644 --- a/c++/descriptor.hpp +++ b/c++/descriptor.hpp @@ -1,6 +1,7 @@ #pragma once #include +#include #include namespace kel { diff --git a/c++/macroscopic.hpp b/c++/macroscopic.hpp index 8a25248..51368e9 100644 --- a/c++/macroscopic.hpp +++ b/c++/macroscopic.hpp @@ -59,5 +59,34 @@ void compute_rho_u ( vel().at({i}) = vel().at({i}) / rho(); } } + +/** + * Calculate the macroscopic variables rho and u in Lattice Units. + */ +template +void compute_rho_u ( + const saw::data>& dfs, + saw::ref> rho, + saw::ref>> vel + ) +{ + using dfi = df_info; + + rho().set(0); + for(uint64_t i = 0; i < Desc::D; ++i){ + vel().at({{i}}).set(0); + } + + for(size_t j = 0; j < Desc::Q; ++j){ + rho() = rho() + dfs(j); + for(size_t i = 0; i < Desc::D; ++i){ + vel().at({{i}}) = vel().at({{i}}) + saw::data{static_cast::type>(dfi::directions[j][i])} * dfs(j); + } + } + + for(size_t i = 0; i < Desc::D; ++i){ + vel().at({{i}}) = vel().at({{i}}) / rho(); + } +} } } diff --git a/c++/util.hpp b/c++/util.hpp index 441ffc8..c5b6a05 100644 --- a/c++/util.hpp +++ b/c++/util.hpp @@ -4,6 +4,8 @@ #include #include +#include +#include namespace kel { namespace lbm { @@ -46,18 +48,37 @@ public: }; */ -void print_progress_bar(uint32_t progress, uint32_t progress_target){ - std::cout - <<"\r" - <<"Progress: " - <<((100 * progress) / progress_target) - <<"% ["; +void print_progress_bar(const uint32_t progress, const uint32_t progress_target){ + std::cout<<"\r"; + // <<"Progress: " + // <<((100 * progress) / progress_target) + // <<"% ["; - uint64_t i{0u}; + uint64_t perc_prog = (100ul*progress) / progress_target; + constexpr uint64_t max_perc_progress = 100u; + perc_prog = std::min(perc_prog, max_perc_progress); - uint64_t max_display = 64u; + std::string_view prog_str{"Progress: "}; + // Progress string + (100 width and perc char) + ([] brackets) + space + pointer + uint64_t i{prog_str.size() + 4u + 2u + 1u + 1u}; + + std::cout< uint64_t{ + struct winsize w; + if(ioctl(STDOUT_FILENO, TIOCGWINSZ,&w) == -1){ + // Standardized Terminal size + return 80u; + } + return w.ws_col; + }(); + max_display = std::max(max_display,i) - i; uint64_t progress_display = (max_display * progress) / progress_target; - for(; i < progress_display; ++i){ + + for(i = 0u; i < progress_display; ++i){ std::cout<<"#"; } for(; i < max_display; ++i){ @@ -65,6 +86,7 @@ void print_progress_bar(uint32_t progress, uint32_t progress_target){ } std::cout<<"]"; + std::cout< #include +#include #include "descriptor.hpp" @@ -55,6 +56,31 @@ struct lbm_vtk_writer> { } }; +template +struct lbm_vtk_writer> { + static saw::error_or apply_header(std::ostream& vtk_file, std::string_view name){ + vtk_file<<"VECTORS "< apply(std::ostream& vtk_file, const saw::data>& field){ + static_assert(D > 0, "Non-dimensionality is bad for velocity."); + static_assert(D <= 3, "4th dimension as well. Mostly due to vtk."); + + // vtk_file<<"VECTORS "< 0){ + vtk_file<<" "; + } + vtk_file< struct lbm_vtk_writer...> , Dim>> { diff --git a/examples/poiseulle_particles_channel_2d.cpp b/examples/poiseulle_particles_channel_2d.cpp index 9f7475d..4673a2a 100644 --- a/examples/poiseulle_particles_channel_2d.cpp +++ b/examples/poiseulle_particles_channel_2d.cpp @@ -54,10 +54,11 @@ using CellStruct = Struct< template using MacroStruct = Struct< - Member, "velocity">, + Member, "velocity">, Member, Member, Member, "force"> + //Member, "debug_acceleration"> >; template @@ -348,13 +349,16 @@ void couple_particles_to_lattice( // p_pos - p_pos_old // auto p_vel_rel = closest_u - p_vel; - saw::data> p_vel_rel; - p_vel_rel.at({{0u}}) = closest_u.at({0u}) - p_vel.at({{0u}}); - p_vel_rel.at({{1u}}) = closest_u.at({1u}) - p_vel.at({{1u}}); + saw::data> p_vel_rel = closest_u - p_vel; + // Add relative velocity to particle acceleration (Since our timestep is 1 we don't need to do anything else) + p_acc = p_acc + p_vel_rel; + + /* for(saw::data i{0u}; i.get() < 2u; ++i){ p_acc.at({{i}}) = p_acc.at({{i}}) + p_vel_rel.at({{i}}); } + */ // Add forces to put away from walls /// 1. Check if neighbour is wall @@ -366,10 +370,11 @@ void couple_particles_to_lattice( auto& n_cell = latt(n_p_cell_pos); auto& n_info = n_cell.template get<"info">()({0u}); - auto& n_macro_cell_particle = macros.at(n_p_cell_pos).template get<"particle">(); + auto& n_macro_cell = macros.at(n_p_cell_pos); + auto& n_macro_cell_particle = n_macro_cell.template get<"particle">(); // If neighbour is wall, then add force pushing the particle away - if(n_info.get() <= 1u or (n_macro_cell_particle.get() != i.get() and n_macro_cell_particle.get() > 0u) ) { + if(n_info.get() <= 1u or (n_macro_cell_particle.get() != (i.get()+1u) and n_macro_cell_particle.get() > 0u) ) { // add to p_acc saw::data dist_dot{ @@ -379,8 +384,16 @@ void couple_particles_to_lattice( if(dist_dot.get() > 0){ // TODO add if particle is close - p_acc.at({{0u}}) = p_acc.at({{0u}}) + saw::data{100 * dfi::directions[dfi::opposite_index[k.get()]][0u]}.template cast_to(); - p_acc.at({{1u}}) = p_acc.at({{1u}}) + saw::data{100 * dfi::directions[dfi::opposite_index[k.get()]][1u]}.template cast_to(); + //auto& n_debug_acc = n_macro_cell.template get<"debug_acceleration">(); + + saw::data> push_force; + { + push_force.at({{0u}}) = saw::data{1000 * dfi::directions[dfi::opposite_index[k.get()]][0u]}.template cast_to(); + push_force.at({{1u}}) = saw::data{1000 * dfi::directions[dfi::opposite_index[k.get()]][1u]}.template cast_to(); + } + + //n_debug_acc = n_debug_acc + push_force; + p_acc = p_acc + push_force; } } } -- cgit v1.2.3