From 34bcca5a8d9234febbed9c11c2c801c3e5d63581 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Thu, 23 Jul 2026 23:20:05 +0200 Subject: Fixed HLBM galilean issue --- modules/core/c++/hlbm.hpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) (limited to 'modules') diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp index 4b1f2b4..bb21ae5 100644 --- a/modules/core/c++/hlbm.hpp +++ b/modules/core/c++/hlbm.hpp @@ -123,6 +123,7 @@ public: auto& part_spheroid_group = part_group; auto& mvel_f = macros.template get<"velocity">(); + auto& mrho_f = macros.template get<"density">(); auto& mpor_f = macros.template get<"porosity">(); auto& particle_N_f = field.template get<"particle_N">(); @@ -160,6 +161,7 @@ public: for(uint64_t i{0u}; i < Desc::D; ++i){ force_p.at({{i}}) = 0.0; } + auto vel_p_old = (pirb_pos-pirb_pos_old) / ts; iterator::apply([&](const auto& index_f) -> void{ // ask for the d_k value here. @@ -193,19 +195,23 @@ public: 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}); + 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& rho = mrho_f.at(index_f); force_p = force_p + momentum; },start,stop); auto& pirb_acc = pirb.template get<"acceleration">(); - pirb_acc = force_p / part_spheroid_group.template get<"total_mass">().at({}); + pirb_acc = - force_p / part_spheroid_group.template get<"total_mass">().at({}); for(saw::data i{0u}; i < sub_steps; ++i){ verlet_step_lambda(pi,ts); -- cgit v1.2.3