From 80deace4d4b2c93945feeec42a8269f3ba2443d2 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Tue, 4 Aug 2026 19:52:08 +0200 Subject: Asking myself how I had a stable implementation before --- modules/core/c++/psm.hpp | 61 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) (limited to 'modules/core/c++/psm.hpp') diff --git a/modules/core/c++/psm.hpp b/modules/core/c++/psm.hpp index 02db1e1..c53dbf9 100644 --- a/modules/core/c++/psm.hpp +++ b/modules/core/c++/psm.hpp @@ -86,8 +86,18 @@ public: template class component final { private: + saw::data relaxation_; + saw::data frequency_; public: - component() = default; + component( + typename saw::native_data_type::type relaxation__ + ): + relaxation_{relaxation__} + { + saw::data one; + one = 1.0; + frequency_ = one / relaxation_; + } template void apply(const saw::data& field, const saw::data& macros, const saw::data& particles, saw::data> index, saw::data time_step) const { @@ -101,6 +111,55 @@ public: auto& rho_f = macros.template get<"density">(); auto& vel_f = macros.template get<"velocity">(); + saw::data>& rho = rho_f.at(index); + saw::data>& vel = vel_f.at(index); + + compute_rho_u(dfs_old_f.at(index),rho,vel); + + auto eq = equilibrium(rho,vel); + + saw::data one{1.0}; + auto& porous = porous_f.at(index); + auto flip_porous = one - porous.at({}); + + auto& dfs = dfs_old_f.at(index); + + auto dfs_cpy = dfs; + + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + uint64_t i_opp = dfi::opposite_index[i]; + dfs.at({i}) = dfs_cpy.at({i}) + frequency_ * (eq.at(i) - dfs_cpy.at({i})) * porous.at({}) + (dfs_cpy.at({i_opp}) - dfs_cpy.at({i}) ) * flip_porous; + } + + auto& force_f = macros.template get<"force">(); + auto& force = force_f.at(index); + for(uint64_t k{0u}; k < Descriptor::D; ++k){ + force.at({{k}}).set(0); + } + + ///////// + saw::data> momentum; + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + + saw::data> e_i; + saw::data> n_ind_i; + for(uint64_t k{0u}; k < Descriptor::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; + } + + // Set Force + force = momentum * porous_f.at(index); } }; -- cgit v1.2.3 From a0767a86817055744230239284441fdbb024fab5 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Wed, 5 Aug 2026 17:13:19 +0200 Subject: WOrking on psm 2way --- modules/core/c++/psm.hpp | 172 +++++++++++++++++++++++++++++++---------------- 1 file changed, 113 insertions(+), 59 deletions(-) (limited to 'modules/core/c++/psm.hpp') diff --git a/modules/core/c++/psm.hpp b/modules/core/c++/psm.hpp index c53dbf9..15ad0f8 100644 --- a/modules/core/c++/psm.hpp +++ b/modules/core/c++/psm.hpp @@ -7,15 +7,30 @@ namespace kel { namespace lbm { namespace cmpt { -struct PSM {}; +struct PsmReset {}; struct PsmOneParticle {}; +struct Psm {}; } +template +class component final { +public: + component() = default; + + template + void apply(const saw::data& field, const saw::data& macros, saw::data> index, saw::data time_step) const { + auto& porosity_f = macros.template get<"porosity">(); + + auto& por = porosity_f.at(index); + por.at({}) = 1.0; + } +}; + /** * PSM collision operator for LBM */ template -class component { +class component { private: saw::data relaxation_; saw::data frequency_; @@ -86,80 +101,119 @@ public: template class component final { private: - saw::data relaxation_; - saw::data frequency_; public: - component( - typename saw::native_data_type::type relaxation__ - ): - relaxation_{relaxation__} - { - saw::data one; - one = 1.0; - frequency_ = one / relaxation_; - } + 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 { + void apply(const saw::data& field, const saw::data& macros, const saw::data& pg, saw::data> index, saw::data time_step, saw::data sub_steps) const { using dfi = df_info; bool is_even = ((time_step.get() % 2) == 0); + + saw::data> one; + one.at({}) = 1.0; auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); - auto& porous_f = macros.template get<"porosity">(); + auto& por_f = macros.template get<"porosity">(); auto& rho_f = macros.template get<"density">(); auto& vel_f = macros.template get<"velocity">(); - - saw::data>& rho = rho_f.at(index); - saw::data>& vel = vel_f.at(index); - - compute_rho_u(dfs_old_f.at(index),rho,vel); - - auto eq = equilibrium(rho,vel); - - saw::data one{1.0}; - auto& porous = porous_f.at(index); - auto flip_porous = one - porous.at({}); - - auto& dfs = dfs_old_f.at(index); - - auto dfs_cpy = dfs; - - for(uint64_t i = 0u; i < Descriptor::Q; ++i){ - uint64_t i_opp = dfi::opposite_index[i]; - dfs.at({i}) = dfs_cpy.at({i}) + frequency_ * (eq.at(i) - dfs_cpy.at({i})) * porous.at({}) + (dfs_cpy.at({i_opp}) - dfs_cpy.at({i}) ) * flip_porous; - } - auto& force_f = macros.template get<"force">(); - auto& force = force_f.at(index); - for(uint64_t k{0u}; k < Descriptor::D; ++k){ - force.at({{k}}).set(0); - } - - ///////// - saw::data> momentum; - for(uint64_t i = 0u; i < Descriptor::Q; ++i){ - - saw::data> e_i; - saw::data> n_ind_i; - for(uint64_t k{0u}; k < Descriptor::D; ++k){ - e_i.at({{k}}) = dfi::directions[i][k]; - n_ind_i.at({k}) = (dfi::directions[i])[k]; + + { + auto parts = pg.template get<"particles">(); + auto parts_size = parts.meta().at({0u}); + + auto& p_coll = pg.template get<"collision">().at({}); + auto& p_rad = p_coll.template get<"radius">(); + + auto& pi = parts.at(index); + auto& pirb = pi.template get<"rigid_body">(); + saw::data>& pirb_pos = pirb.template get<"position">(); + auto& pirb_pos_old = pirb.template get<"position_old">(); + + saw::data> ts; + ts.at({}) = one.at({}) / sub_steps.template cast_to(); + + saw::data> start; + saw::data> stop; + + auto eo_aabb = particle_aabb::calculate(pg,{{0u}},vel_f.meta()); + if(eo_aabb.is_error()){ + return; } + auto& aabb = eo_aabb.get_value(); - uint64_t i_opp = dfi::opposite_index[i]; + /// Ok, I iterate over the space which covers our particle? So lower bounds to upper bounds + start = aabb.template get<"a">(); + stop = aabb.template get<"b">(); - 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; + saw::data> force_p; + for(uint64_t i{0u}; i < Descriptor::D; ++i){ + force_p.at({{i}}) = 0.0; + } + auto vel_p_old = (pirb_pos-pirb_pos_old) / ts; - momentum = momentum + ei_dfs; + iterator::apply([&](const auto& index_f) -> void{ + // ask for the d_k value here. + // For every value im iterating over I need sth + // std::cout<<"Pos: "<> rel_dist = saw::math::vectorize_data(index_f).template cast_to() - pirb_pos; + saw::data> eps; + eps.at({}) = 1.5f; + + auto& por = por_f.at(index_f); + por = particle_porosity>::calculate(rel_dist,p_rad,eps); + + if(por.at({}).get() >= 1.0f){ + return; + } + + saw::data> momentum; + for(uint64_t i{0u}; i < Descriptor::D; ++i){ + momentum.at({{i}}) = 0.0; + } + + for(uint64_t i{0u}; i < Descriptor::Q; ++i){ + saw::data> e_i; + saw::data> n_ind_i; + for(uint64_t k{0u}; k < Descriptor::D; ++k){ + e_i.at({{k}}) = (dfi::directions[i])[k]; + n_ind_i.at({k}) = index_f.at({k}) + (dfi::directions[i])[k]; + } + + uint64_t i_opp = dfi::opposite_index[i]; + + auto u_p_e = saw::math::dot(e_i,vel_p_old); + + saw::data dfs_added = dfs.at({i})*(saw::data{1}-u_p_e.at({})) + dfs_old_f.at(n_ind_i).at({i_opp})*(saw::data{1}+u_p_e.at({})); + saw::data> dfs_added_v; + dfs_added_v.at({}) = dfs_added; + auto ei_dfs = e_i * dfs_added_v; + + momentum = momentum + ei_dfs; + } + // technically needs to adjust for rotation as well + + auto& force = force_f.at(index_f); + auto& rho = rho_f.at(index_f); + // To Fluid + + auto flip_por = one - por; + force = momentum * flip_por; + // To Particle + force_p = force_p - force; + },start,stop); + + auto& pirb_acc = pirb.template get<"acceleration">(); + pirb_acc = force_p / pg.template get<"total_mass">().at({}); + + for(saw::data i{0u}; i < sub_steps; ++i){ + verlet_step_lambda(pi,ts); + } } - - // Set Force - force = momentum * porous_f.at(index); } }; -- cgit v1.2.3