From 1693f8ccb3fbdbc8a1bd35a5ffde684ca4c68a3f Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Sat, 18 Jul 2026 21:45:14 +0200 Subject: Dangling --- .../moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp | 29 ++++++++++- modules/core/c++/hlbm.hpp | 59 +++++++++++++++++----- modules/core/c++/particle/particle.hpp | 39 +++++++++++++- modules/core/c++/particle/porosity.hpp | 6 +-- modules/core/c++/psm.hpp | 25 +++++++++ 5 files changed, 140 insertions(+), 18 deletions(-) diff --git a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp index a4cae03..a41d2b7 100644 --- a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp +++ b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp @@ -146,7 +146,7 @@ saw::error_or setup_initial_conditions( ); saw::data> rad; - rad.at({}) = 2.0f; + rad.at({}) = dim_y * 0.1f; saw::data> dense; dense.at({}) = 1.0f; particles = create_spheroid_particle_group(rad,dense,{16u}); @@ -173,6 +173,18 @@ saw::error_or step( auto& info_f = fields.template get<"info">(); auto& porous_f = macros.template get<"porosity">(); + q.submit([&](acpp::sycl::handler& h){ + component> hlbm_reset; + h.parallel_for(acpp::sycl::range{dim_x,dim_y}, [=](acpp::sycl::id idx){ + saw::data> index; + for(uint64_t i = 0u; i < Desc::D; ++i){ + index.at({{i}}).set(idx[i]); + } + + hlbm_reset.apply(fields,macros,index,t_i); + }); + }).wait(); + q.submit([&](acpp::sycl::handler& h){ component> hlbm_one_part; @@ -361,6 +373,21 @@ saw::error_or lbm_main(int argc, char** argv){ } sycl_q.wait(); + if( i.get() % 128u == 0u){ + { + auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_vtk_file(out_dir,"m",i.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + } + // Stream sycl_q.submit([&](acpp::sycl::handler& h){ component> stream; diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp index cc35f9b..25b6bf0 100644 --- a/modules/core/c++/hlbm.hpp +++ b/modules/core/c++/hlbm.hpp @@ -11,7 +11,7 @@ namespace kel { namespace lbm { namespace cmpt { -struct HlbmInit {}; +struct HlbmReset {}; struct Hlbm {}; struct HlbmOneParticle {}; struct HlbmParticle {}; @@ -19,15 +19,9 @@ struct HlbmOneParticleMomentumExchange {}; } template -class component final { -private: - typename saw::native_data_type::type relaxation_; - saw::data frequency_; +class component final { public: - component(typename saw::native_data_type::type relaxation__): - relaxation_{relaxation__}, - frequency_{typename saw::native_data_type::type(1) / relaxation_} - {} + component() = default; template void apply(const saw::data& field, const saw::data& macros, saw::data> index, saw::data time_step) const { @@ -36,13 +30,13 @@ public: auto& particle_D_f = field.template get<"particle_D">(); auto& por = porosity_f.at(index); - por = {}; + por.at({}) = 1.0; auto& pnf = particle_N_f.at(index); pnf = {}; auto& pnd = particle_D_f.at(index); - pnd = {}; + pnd.at({}) = 0.0; } }; @@ -96,8 +90,6 @@ public: dfs_old_f.at(index).at({i}) = dfs_old_f.at(index).at({i}) + frequency_ * (eq.at(i) - dfs_old_f.at(index).at({i})); } - D.at({}) = 0.0; - N = {}; } }; @@ -118,13 +110,20 @@ public: /// Iterate over the grid bounds // auto& grid = p.template get<"grid">(); + using dfi = df_info; + bool is_even = ((time_step.get() % 2) == 0); + + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); auto& part_spheroid_group = part_group; auto& mvel = macros.template get<"velocity">(); + auto& mpor_f = macros.template get<"porosity">(); { auto& parts = part_spheroid_group.template get<"particles">(); auto parts_size = parts.meta().at({0u}); + auto& p_coll = part_spheroid_group.template get<"collision">().at({}); + auto& p_rad = p_coll.template get<"radius">(); auto& pi = parts.at(index); auto& pirb = pi.template get<"rigid_body">(); @@ -143,14 +142,48 @@ public: start = aabb.template get<"a">(); stop = aabb.template get<"b">(); + saw::data> force_p; + iterator::apply([&](const auto& index){ // ask for the d_k value here. // For every value im iterating over I need sth // std::cout<<"Pos: "<> momentum; + for(uint64_t i = 0u; i < Desc::Q; ++i){ + saw::data> e_i; + saw::data> n_ind_i; + for(uint64_t k{0u}; k < Desc::D; ++k){ + e_i.at({{k}}) = (dfi::directions[i])[k]; + n_ind_i.at({k}) = (dfi::directions[i])[k]; + } + + uint64_t i_opp = dfi::opposite_index[i]; + + saw::data dfs_added = dfs.at({i}) - dfs_old_f.at(n_ind_i).at({i_opp}); + saw::data> dfs_added_v; + dfs_added_v.at({}) = dfs_added; + auto ei_dfs = e_i * dfs_added_v; + + momentum = momentum + ei_dfs; + } + + auto& mpor = mpor_f.at(index); + auto rel_dist = saw::math::vectorize_data(index).template cast_to() - pirb_pos; + saw::data> eps; + eps.at({}) = 1.5f; + mpor = particle_porosity>::calculate(rel_dist,p_rad,eps); + force_p = force_p + momentum * mpor; },start,stop); + auto& pirb_acc = pirb.template get<"acceleration">(); + pirb_acc = force_p; + + verlet_step_lambda(part_spheroid_group,{{0u}},{{1u}}); + // Check } } diff --git a/modules/core/c++/particle/particle.hpp b/modules/core/c++/particle/particle.hpp index 3c3630a..6d57930 100644 --- a/modules/core/c++/particle/particle.hpp +++ b/modules/core/c++/particle/particle.hpp @@ -137,7 +137,44 @@ saw::data>> create_spheroid */ template -constexpr auto verlet_step_lambda = [](saw::data>& particle, saw::data> time_step_delta){ +constexpr auto verlet_step_lambda_old = [](saw::data>& particle, saw::data> time_step_delta){ + auto& body = particle.template get<"rigid_body">(); + + auto& pos = body.template get<"position">(); + auto& pos_old = body.template get<"position_old">(); + + auto& pos_acc = body.template get<"acceleration">(); + + auto& rot = body.template get<"rotation">(); + auto& rot_old = body.template get<"rotation_old">(); + + auto& rot_acc = body.template get<"angular_acceleration">(); + + auto tsd_squared = time_step_delta * time_step_delta; + + saw::data> pos_new; + // Actual step + saw::data> two; + two.at({}).set(2.0); + pos_new = pos * two - pos_old + pos_acc * tsd_squared; + + // Angular + saw::data::Schema> rot_new; + rot_new = rot * two - rot_old + rot_acc * tsd_squared; + + // Swap - Could be std::swap? + pos_old = pos; + pos = pos_new; + + rot_old = rot; + rot = rot_new; +}; + +template +constexpr auto verlet_step_lambda = [](saw::data>& pgrp, saw::data> index, saw::data> time_step_delta){ + auto& parts = pgrp.template get<"particles">(); + auto& particle = parts.at(index); + auto& body = particle.template get<"rigid_body">(); auto& pos = body.template get<"position">(); diff --git a/modules/core/c++/particle/porosity.hpp b/modules/core/c++/particle/porosity.hpp index 353cd58..4d95cce 100644 --- a/modules/core/c++/particle/porosity.hpp +++ b/modules/core/c++/particle/porosity.hpp @@ -86,9 +86,9 @@ public: saw::data> dist_len = saw::math::dot(dist,dist); dist_len.at({}).set(std::sqrt(dist_len.at({}).get())); - auto& coll = part_group.template get<"collision">(); + auto& coll = part_group.template get<"collision">().at({}); const saw::data>& rad_d = coll.template get<"radius">(); - const auto& eps = part_group.template get<"epsilon">(); + const auto& eps = part_group.template get<"epsilon">().at({}); // Move this somewhere saw::data> eps_h; @@ -116,7 +116,7 @@ public: * cos^2 ( ||x-X(t)||_2 - (R-eps/2) ) */ { - typename saw::native_data_type::type inner = (std::numbers::pi / (eps)) * (dist_len - rad_d_eps_n).at({}).get(); + typename saw::native_data_type::type inner = (std::numbers::pi / (eps.at({}).get())) * (dist_len - rad_d_eps_n).at({}).get(); auto cos_inner = std::cos(inner); auto cos_i_2 = cos_inner * cos_inner; inside.at({}).set(cos_i_2); diff --git a/modules/core/c++/psm.hpp b/modules/core/c++/psm.hpp index 30a00da..02db1e1 100644 --- a/modules/core/c++/psm.hpp +++ b/modules/core/c++/psm.hpp @@ -8,6 +8,7 @@ namespace kel { namespace lbm { namespace cmpt { struct PSM {}; +struct PsmOneParticle {}; } /** @@ -79,5 +80,29 @@ public: } }; +/** + * PSM collision operator for LBM + */ +template +class component final { +private: +public: + component() = default; + + template + void apply(const saw::data& field, const saw::data& macros, const saw::data& particles, saw::data> index, saw::data time_step) const { + + using dfi = df_info; + bool is_even = ((time_step.get() % 2) == 0); + + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + auto& porous_f = macros.template get<"porosity">(); + + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"velocity">(); + + } +}; + } } -- cgit v1.2.3