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++/fplbm.hpp | 126 +++++++++++++++++++++++++++++++++++++++------ 1 file changed, 110 insertions(+), 16 deletions(-) (limited to 'modules/core/c++/fplbm.hpp') diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp index cfb0062..4acdd36 100644 --- a/modules/core/c++/fplbm.hpp +++ b/modules/core/c++/fplbm.hpp @@ -10,6 +10,7 @@ namespace cmpt { struct FpLbmReset{}; struct FpLbm {}; struct FpLbmOneParticle {}; +struct FpLbmOneParticleImplicit {}; struct FpLbmOneParticleNoVelocity{}; } namespace method { @@ -36,18 +37,20 @@ public: for(uint64_t i{0u}; i < Descriptor::D; ++i){ force.at({{i}}) = 0.0; } + auto info = field.template get<"info">().at(index).get(); + if(info == 2u){ + bool is_even = ((time_step.get() % 2) == 0); - bool is_even = ((time_step.get() % 2) == 0); - - auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); - auto& dfs = dfs_old_f.at(index); + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + auto& dfs = dfs_old_f.at(index); - auto& rho_f = macros.template get<"density">(); - auto& rho = rho_f.at(index); - auto& vel_f = macros.template get<"velocity">(); - auto& vel = vel_f.at(index); + auto& rho_f = macros.template get<"density">(); + auto& rho = rho_f.at(index); + auto& vel_f = macros.template get<"velocity">(); + auto& vel = vel_f.at(index); - compute_rho_u(dfs,rho,vel); + compute_rho_u(dfs,rho,vel); + } } }; @@ -83,7 +86,7 @@ public: auto& vel_f = macros.template get<"velocity">(); saw::data>& vel = vel_f.at(index); - compute_rho_u(dfs,rho,vel); + // compute_rho_u(dfs,rho,vel); auto eq = equilibrium(rho,vel); using dfi = df_info; @@ -158,10 +161,94 @@ public: auto& pirb_pos = pirb.template get<"position">(); auto& pirb_pos_old = pirb.template get<"position_old">(); + + auto& p_coll = pg.template get<"collision">().at({}); + auto& p_rad = p_coll.template get<"radius">(); + + auto eo_aabb = particle_aabb::calculate(pg,{{0u}},vel_f.meta()); + if(eo_aabb.is_error()){ + return; + } + auto& aabb = eo_aabb.get_value(); + + saw::data> two; + two.at({}).set(2.0f); + + saw::data> one; + one.at({}) = 1.0; + + saw::data> eps; + eps.at({}) = 1.5f; + saw::data> sss; - sss.at({}) = sub_steps.template cast_to(); + sss.at({}) = one.at({}) / sub_steps.template cast_to(); auto vel_s = (pirb_pos-pirb_pos_old) / sss; + saw::data> force_p{}; + + iterator::apply([&](const auto& index_f) -> void { + auto& force = force_f.at(index_f); + auto& vel = vel_f.at(index_f); + auto& por = por_f.at(index_f); + auto& rho = rho_f.at(index_f); + + saw::data> rel_dist = saw::math::vectorize_data(index_f).template cast_to() - pirb_pos; + + por = particle_porosity>::calculate(rel_dist,p_rad,eps); + auto flip_por = one - por; + + // vel_s is technically time the density of the particle? + + force = ( vel_s * rho - vel * rho ) * two * flip_por; + + force_p = force_p - force; + }, aabb.template get<"a">(), aabb.template get<"b">()); + + 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,sss); + } + } +}; + +template +class component final { +public: + using Component = cmpt::FpLbmOneParticle; +private: +public: + component() = default; + + template + 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 { + // void apply(saw::data& field, saw::data> index, saw::data time_step){ + bool is_even = ((time_step.get() % 2) == 0); + + //auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + // auto& dfs = dfs_old_f.at(index); + + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"velocity">(); + auto& por_f = macros.template get<"porosity">(); + auto& force_f = macros.template get<"force">(); + + /** + * The other methods all work with a flipped porosity. + * For compat reasons this is also flipped + */ + // TODO - Change to tuple later + auto parts = pg.template get<"particles">(); + + auto parts_size = parts.meta().at({0u}); + auto& pi = parts.at(index); + + auto& pirb = pi.template get<"rigid_body">(); + auto& pirb_pos = pirb.template get<"position">(); + auto& pirb_pos_old = pirb.template get<"position_old">(); + + auto& p_coll = pg.template get<"collision">().at({}); auto& p_rad = p_coll.template get<"radius">(); @@ -172,15 +259,19 @@ public: auto& aabb = eo_aabb.get_value(); saw::data> two; - two.at({}).set(2); + two.at({}).set(2.0f); saw::data> one; one.at({}) = 1.0; saw::data> eps; eps.at({}) = 1.5f; + + saw::data> sss; + sss.at({}) = one.at({}) / sub_steps.template cast_to(); + auto vel_s = (pirb_pos-pirb_pos_old) / sss; - saw::data> force_p; + saw::data> force_p{}; iterator::apply([&](const auto& index_f) -> void { auto& force = force_f.at(index_f); @@ -193,13 +284,16 @@ public: por = particle_porosity>::calculate(rel_dist,p_rad,eps); auto flip_por = one - por; - force = ( vel_s * rho - vel * rho ) * two * flip_por; + // vel_s is technically time the density of the particle? + + force = ( vel_s * rho - vel * rho ) * flip_por / (one + flip_por / two); + // force = ( vel_s * rho - vel * rho ) * two * flip_por; - force_p = force_p + force; + force_p = force_p - force; }, aabb.template get<"a">(), aabb.template get<"b">()); auto& pirb_acc = pirb.template get<"acceleration">(); - pirb_acc = - force_p / pg.template get<"total_mass">().at({}); + pirb_acc = force_p / pg.template get<"total_mass">().at({}); for(saw::data i{0u}; i < sub_steps; ++i){ verlet_step_lambda(pi,sss); -- 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++/fplbm.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) (limited to 'modules/core/c++/fplbm.hpp') diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp index 4acdd36..1f38f51 100644 --- a/modules/core/c++/fplbm.hpp +++ b/modules/core/c++/fplbm.hpp @@ -208,7 +208,7 @@ public: pirb_acc = force_p / pg.template get<"total_mass">().at({}); for(saw::data i{0u}; i < sub_steps; ++i){ - //verlet_step_lambda(pi,sss); + verlet_step_lambda(pi,sss); } } }; @@ -286,8 +286,7 @@ public: // vel_s is technically time the density of the particle? - force = ( vel_s * rho - vel * rho ) * flip_por / (one + flip_por / two); - // force = ( vel_s * rho - vel * rho ) * two * flip_por; + force = ( vel_s * rho - vel * rho ) * two * flip_por / (one + flip_por); force_p = force_p - force; }, aabb.template get<"a">(), aabb.template get<"b">()); -- cgit v1.2.3