From 51a044e3160df05ad56102d3b8b1e0087c60d111 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Mon, 4 Aug 2025 16:16:17 +0200 Subject: emits weird behaviour on particles --- examples/poiseulle_particles_channel_2d.cpp | 65 ++++++++++++++++++++--------- 1 file changed, 46 insertions(+), 19 deletions(-) (limited to 'examples') diff --git a/examples/poiseulle_particles_channel_2d.cpp b/examples/poiseulle_particles_channel_2d.cpp index 9003d1b..725e008 100644 --- a/examples/poiseulle_particles_channel_2d.cpp +++ b/examples/poiseulle_particles_channel_2d.cpp @@ -6,6 +6,8 @@ #include +#include + namespace kel { namespace lbm { @@ -45,7 +47,7 @@ using CellStruct = Struct< Member, "dfs">, Member, "dfs_old">, Member, "info">, - Member, "particle_mask">, + // Member, "particle_mask">, Member, "forcing"> >; @@ -293,7 +295,7 @@ void add_particles(kel::lbm::particle_system& part_sys){ auto& p_size = part.template get<"size">(); { auto& pos = rigid_body.template get<"position">(); - auto& old_pos = rigid_body.template get<"position">(); + auto& old_pos = rigid_body.template get<"position_old">(); pos = {{32u,64u}}; old_pos = {{32u,64u}}; @@ -302,7 +304,7 @@ void add_particles(kel::lbm::particle_system& part_sys){ } { - auto eov = particle_sys.add(part); + auto eov = part_sys.add_particle(part); if(eov.is_error()){ exit(-1); } @@ -310,9 +312,9 @@ void add_particles(kel::lbm::particle_system& part_sys){ } void couple_particles_to_lattice( - kel::lbm::particle_system& part_sys, + kel::lbm::particle_system& part_sys, saw::data& latt, - saw::data,kel::lbm::sch::D2Q9::D>>& macros, + saw::data,kel::lbm::sch::D2Q9::D>>& macros, saw::data time_step ){ using namespace kel::lbm; @@ -322,7 +324,7 @@ void couple_particles_to_lattice( iterate_over([&](const saw::data>& index){ auto& cell = latt(index); auto& info = cell.template get<"info">(); - auto& mask = cell.template get<"particle_mask">(); + // auto& mask = cell.template get<"particle_mask">(); for(saw::data i{0u}; i < part_sys.size(); ++i){ auto& part = part_sys.at(i); @@ -350,17 +352,21 @@ void couple_particles_to_lattice( // Spread force to close fluid cells iterate_over([&](const saw::data>& index){ (void)index; - }, {{0u,0u}}, p_mask.dims()); + }, {{0u,0u}}, p_mask.template get<"grid">().dims()); // Fluid to Particle Coupling iterate_over([&](const saw::data>& index){ saw::data> mask_shift{{ - index.at({0u}) * x_dir.at({0u}) + index.at({0u}) * y_dir.at({0u}), - index.at({1u}) * x_dir.at({1u}) + index.at({1u}) * y_dir.at({1u}) + index.at({0u}).template cast_to() * x_dir.at({0u}) + index.at({0u}).template cast_to() * y_dir.at({0u}), + index.at({1u}).template cast_to() * x_dir.at({1u}) + index.at({1u}).template cast_to() * y_dir.at({1u}) }}; // Lagrange in Euler - auto p_pos_lie = p_pos + mask_shift; + saw::data> p_pos_lie {{ + p_pos.at({0u}) + mask_shift.at({0u}), + p_pos.at({1u}) + mask_shift.at({1u}) + }}; + // Cast down to get lower corner. // Before casting shift by 0.5 for closest pick saw::data> p_cell_pos {{ @@ -371,19 +377,17 @@ void couple_particles_to_lattice( // For now pick the closest U auto& closest_u = macros.at(p_cell_pos).template get<"velocity">(); auto p_shift = closest_u; - for(saw::data i{0u}; i < 2u; ++i){ - p_shift.at(i) = p_shift.at(i) * time_step; + for(saw::data i{0u}; i.get() < 2u; ++i){ + p_pos.at(i) = p_pos.at(i) + p_shift.at(i) * time_step; } - - // Add forces to put away from walls /// 1. Check if neighbour is wall /// 2. Add force pushing away from wall // Add forces to push away from other particles - }, {{0u,0u}}, p_mask.dims()); + }, {{0u,0u}}, p_mask.template get<"grid">().dims()); } }, {{0u,0u}}, meta); @@ -514,14 +518,14 @@ int main(int argc, char** argv){ set_geometry(lattice); { - saw::data, sch::D2Q9::D>> geo{dim}; iterate_over([&](const saw::data>& index){ + auto& cell = lattice(index); auto& info = cell.template get<"info">(); - geo(index).template get<"info">().set(info({0u}).get()); + }, {{0u,0u}}, dim); write_vtk_file(out_dir / "geometry.vtk", geo); @@ -534,7 +538,7 @@ int main(int argc, char** argv){ saw::data,sch::D2Q9::D>> macros{dim}; - for(uint64_t i = 0u; i < 4096u*16u; ++i){ + for(uint64_t i = 0u; i < 4096u*4u; ++i){ bool even_step = ((i % 2u) == 0u); { @@ -548,10 +552,32 @@ int main(int argc, char** argv){ auto& part_mask = macros.at(index).template get<"particle">(); compute_rho_u(dfs,rho,vel); rho = rho * saw::data{dfi::cs2}; - part = cell.template get<"particle_mask">()({0u}); + + part_mask.set(0u); }, {{0u,0u}}, meta); + for(saw::data i{0u}; i < particle_sys.size(); ++i){ + auto& p = particle_sys.at({i}); + auto& p_rb = p.template get<"rigid_body">(); + auto& p_pos = p_rb.template get<"position">(); + + saw::data> p_cell_pos {{ + static_cast(p_pos.at({0u}).get()+0.5), + static_cast(p_pos.at({1u}).get()+0.5) + }}; + + std::cout<<"Pos: "<= dim.at({0u}) ){ + p_cell_pos.at({0u}).set(dim.at({0u}).get() - 1u ); + } + if(p_cell_pos.at({1u}) >= dim.at({1u}) ){ + p_cell_pos.at({1u}).set(dim.at({1u}).get() - 1u ); + } + + macros(p_cell_pos).template get<"particle">().set(i.get()); + } { std::string vtk_f_name{"macros_"}; @@ -562,6 +588,7 @@ int main(int argc, char** argv){ couple_particles_to_lattice(particle_sys, lattice, macros, {1u}); lbm_step(lattice, i); + particle_sys.step({1u}); } return 0; -- cgit v1.2.3