From 7220084ddae1f4c6007f5111913cb7f87a794657 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Thu, 13 Aug 2026 13:11:05 +0200 Subject: Renaming velocity to momentum. for now --- modules/core/c++/boundary.hpp | 36 +-- modules/core/c++/collision.hpp | 4 +- modules/core/c++/component/psm.hpp | 2 +- modules/core/c++/fplbm.hpp | 8 +- modules/core/c++/hlbm.hpp | 4 +- modules/core/c++/macroscopic.hpp | 1 + modules/core/c++/particle.hpp | 2 +- modules/core/c++/psm.hpp | 516 +++++++++++++++++++++++++++++-------- modules/core/c++/write_csv.hpp | 2 +- modules/core/c++/write_vtk.hpp | 2 +- 10 files changed, 440 insertions(+), 137 deletions(-) (limited to 'modules/core/c++') diff --git a/modules/core/c++/boundary.hpp b/modules/core/c++/boundary.hpp index 7f185ad..e7567c9 100644 --- a/modules/core/c++/boundary.hpp +++ b/modules/core/c++/boundary.hpp @@ -133,14 +133,14 @@ template class component final { private: saw::data> density_; - saw::data> velocity_; + saw::data> momentum_; public: component( saw::data> density__, - saw::data> velocity__ + saw::data> momentum__ ): density_{density__}, - velocity_{velocity__} + momentum_{momentum__} {} template @@ -151,7 +151,7 @@ public: auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); - auto eq = equilibrium(density_,velocity_); + auto eq = equilibrium(density_,momentum_); dfs_old_f.at(index) = eq; } @@ -726,12 +726,12 @@ public: template class component, Encode> final { private: - saw::data> velocity_; + saw::data> momentum_; public: component( - saw::data> velocity__ + saw::data> momentum__ ): - velocity_{velocity__} + momentum_{momentum__} {} template @@ -753,20 +753,20 @@ public: } saw::data> rho; - rho.at({}) = (dfs.at({0u}) + dfs.at({4u}) + dfs.at({3u}) + dir_sum * 2) / (velocity_.at({{0u}}) + saw::data{static_cast::type>(1.0)}); + rho.at({}) = (dfs.at({0u}) + dfs.at({4u}) + dfs.at({3u}) + dir_sum * 2) / (momentum_.at({{0u}}) + saw::data{static_cast::type>(1.0)}); if constexpr (East) { - dfs.at({2u}) = dfs.at({1u}) + saw::data{static_cast::type>(2.0 / 3.0)} * rho.at({}) * velocity_.at({{0u}}); - dfs.at({6u}) = dfs.at({5u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}}) - + (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * velocity_.at({{1u}})) * 0.5f; - dfs.at({8u}) = dfs.at({7u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}}) - + (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * velocity_.at({{1u}})) * 0.5f; + dfs.at({2u}) = dfs.at({1u}) + saw::data{static_cast::type>(2.0 / 3.0)} * rho.at({}) * momentum_.at({{0u}}); + dfs.at({6u}) = dfs.at({5u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}}) + + (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * momentum_.at({{1u}})) * 0.5f; + dfs.at({8u}) = dfs.at({7u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}}) + + (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * momentum_.at({{1u}})) * 0.5f; }else if constexpr ( not East ){ - dfs.at({1u}) = dfs.at({2u}) + saw::data{static_cast::type>(2.0 / 3.0)} * rho.at({}) * velocity_.at({{0u}}); - dfs.at({5u}) = dfs.at({6u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}}) - + (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * velocity_.at({{1u}}))*0.5f; - dfs.at({7u}) = dfs.at({8u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}}) - + (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * velocity_.at({{1u}}))*0.5f; + dfs.at({1u}) = dfs.at({2u}) + saw::data{static_cast::type>(2.0 / 3.0)} * rho.at({}) * momentum_.at({{0u}}); + dfs.at({5u}) = dfs.at({6u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}}) + + (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * momentum_.at({{1u}}))*0.5f; + dfs.at({7u}) = dfs.at({8u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}}) + + (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * momentum_.at({{1u}}))*0.5f; } } }; diff --git a/modules/core/c++/collision.hpp b/modules/core/c++/collision.hpp index 023f61f..39514b8 100644 --- a/modules/core/c++/collision.hpp +++ b/modules/core/c++/collision.hpp @@ -75,7 +75,7 @@ public: auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data>& rho = rho_f.at(index); saw::data>& vel = vel_f.at(index); @@ -135,7 +135,7 @@ public: auto& dfs = dfs_old_f.at(index); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data>& rho = rho_f.at(index); diff --git a/modules/core/c++/component/psm.hpp b/modules/core/c++/component/psm.hpp index cada53c..5832760 100644 --- a/modules/core/c++/component/psm.hpp +++ b/modules/core/c++/component/psm.hpp @@ -39,7 +39,7 @@ public: auto& porous_f = macros.template get<"porosity">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data>& rho = rho_f.at(index); saw::data>& vel = vel_f.at(index); diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp index e2527ce..382d067 100644 --- a/modules/core/c++/fplbm.hpp +++ b/modules/core/c++/fplbm.hpp @@ -46,7 +46,7 @@ public: auto& rho_f = macros.template get<"density">(); auto& rho = rho_f.at(index); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& vel = vel_f.at(index); compute_rho_u(dfs,rho,vel); @@ -86,7 +86,7 @@ public: auto& rho_f = macros.template get<"density">(); saw::data>& rho = rho_f.at(index); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data> vel = vel_f.at(index); auto& force_f = macros.template get<"force">(); @@ -148,7 +148,7 @@ public: // auto& dfs = dfs_old_f.at(index); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); auto& force_f = macros.template get<"force">(); @@ -245,7 +245,7 @@ public: auto& dfs = dfs_old_f.at(index); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); // Temporary find a better way diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp index cc21dcf..e8ea66a 100644 --- a/modules/core/c++/hlbm.hpp +++ b/modules/core/c++/hlbm.hpp @@ -70,7 +70,7 @@ public: auto& porosity_f = macros.template get<"porosity">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data>& rho = rho_f.at(index); saw::data>& vel = vel_f.at(index); @@ -125,7 +125,7 @@ public: auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); auto& part_spheroid_group = part_group; - auto& mvel_f = macros.template get<"velocity">(); + auto& mvel_f = macros.template get<"momentum">(); auto& mrho_f = macros.template get<"density">(); auto& mpor_f = macros.template get<"porosity">(); diff --git a/modules/core/c++/macroscopic.hpp b/modules/core/c++/macroscopic.hpp index 19ab3e9..768fef0 100644 --- a/modules/core/c++/macroscopic.hpp +++ b/modules/core/c++/macroscopic.hpp @@ -28,6 +28,7 @@ void compute_rho_u ( } } + // TODO - Change to momentum, but that needs checks in every impacted function for(size_t i = 0; i < Desc::D; ++i){ vel().at({{i}}) = vel().at({{i}}) / rho().at({}); } diff --git a/modules/core/c++/particle.hpp b/modules/core/c++/particle.hpp index 2c60884..efcc3fb 100644 --- a/modules/core/c++/particle.hpp +++ b/modules/core/c++/particle.hpp @@ -25,7 +25,7 @@ public: // Compute forces - // Update particle velocity + // Update particle momentum verlet_step_lambda_old(p,{1.0}); // Update porosity over lattice nodes diff --git a/modules/core/c++/psm.hpp b/modules/core/c++/psm.hpp index 15ad0f8..2daff5e 100644 --- a/modules/core/c++/psm.hpp +++ b/modules/core/c++/psm.hpp @@ -55,7 +55,7 @@ public: auto& porous_f = macros.template get<"porosity">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data>& rho = rho_f.at(index); saw::data>& vel = vel_f.at(index); @@ -105,116 +105,418 @@ 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 { - - using dfi = df_info; - bool is_even = ((time_step.get() % 2) == 0); - - saw::data> one; - one.at({}) = 1.0; +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 { - auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + using dfi = df_info; - auto& por_f = macros.template get<"porosity">(); - auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); - auto& force_f = macros.template get<"force">(); - - { - 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(); + bool is_even = ((time_step.get() % 2u) == 0u); - /// 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">(); + /* + * ------------------------------------------------------------ + * Fields + * ------------------------------------------------------------ + */ - 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; - - 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); - } - } - } + auto& dfs_old_f = + (is_even) + ? field.template get<"dfs_old">() + : field.template get<"dfs">(); + + auto& por_f = macros.template get<"porosity">(); + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"momentum">(); + auto& force_f = macros.template get<"force">(); + + saw::data> one; + one.at({}) = 1.0; + + saw::data> half; + half.at({}) = 0.5; + + /* + * ------------------------------------------------------------ + * Particle + * ------------------------------------------------------------ + */ + + { + auto parts = pg.template get<"particles">(); + + 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">(); + + + /* + * -------------------------------------------------------- + * Particle timestep + * -------------------------------------------------------- + */ + + saw::data> ts; + + ts.at({}) = + one.at({}) + / sub_steps.template cast_to(); + + + /* + * -------------------------------------------------------- + * Particle momentum + * -------------------------------------------------------- + */ + + auto vel_p_old = + (pirb_pos - pirb_pos_old) / ts; + + + /* + * -------------------------------------------------------- + * Find particle AABB + * -------------------------------------------------------- + */ + + auto eo_aabb = + particle_aabb< + typename ParticleSchema::ValueType + >::calculate( + pg, + {{0u}}, + vel_f.meta() + ); + + if(eo_aabb.is_error()){ + return; + } + + auto& aabb = eo_aabb.get_value(); + + saw::data> start; + saw::data> stop; + + start = aabb.template get<"a">(); + stop = aabb.template get<"b">(); + + + /* + * -------------------------------------------------------- + * Total force acting on particle + * -------------------------------------------------------- + */ + + saw::data> force_p; + + for(uint64_t k = 0u; k < Descriptor::D; ++k){ + force_p.at({{k}}) = 0.0; + } + + + /* + * -------------------------------------------------------- + * Iterate over particle AABB + * -------------------------------------------------------- + */ + + iterator::apply( + [&](const auto& index_f) -> void { + + /* + * ------------------------------------------------ + * Distribution function at current cell + * ------------------------------------------------ + */ + + auto& dfs = + dfs_old_f.at(index_f); + + + /* + * ------------------------------------------------ + * Particle relative position + * ------------------------------------------------ + */ + + saw::data> rel_dist = + saw::math::vectorize_data(index_f) + .template cast_to() + - pirb_pos; + + + /* + * ------------------------------------------------ + * PSM particle porosity + * ------------------------------------------------ + */ + + saw::data> eps; + eps.at({}) = 1.5; + + auto& por = + por_f.at(index_f); + + por = + particle_porosity< + T, + Descriptor::D, + 1u, + por::ParticleSpheroid + >::calculate( + rel_dist, + p_rad, + eps + ); + + + /* + * Outside particle. + */ + + if(por.at({}).get() >= 1.0){ + return; + } + + + /* + * ------------------------------------------------ + * Momentum exchange for this cell + * ------------------------------------------------ + */ + + saw::data> momentum; + + for(uint64_t k = 0u; k < Descriptor::D; ++k){ + momentum.at({{k}}) = 0.0; + } + + + /* + * ------------------------------------------------ + * D3Q27 / general Q momentum exchange + * + * We exclude i=0 because the rest population has + * no momentum. + * + * The original code summed both directions of + * every opposite pair. Therefore we apply 1/2 + * after the sum. + * ------------------------------------------------ + */ + + for(uint64_t i = 1u; i < Descriptor::Q; ++i){ + + uint64_t i_opp = + dfi::opposite_index[i]; + + + /* + * Direction vector e_i. + */ + + saw::data> e_i; + + saw::data< + sch::FixedArray< + sch::UInt64, + Descriptor::D + > + > n_ind_i; + + + for(uint64_t k = 0u; + k < Descriptor::D; + ++k){ + + e_i.at({{k}}) = + dfi::directions[i][k]; + + /* + * Neighbor in direction e_i. + */ + + n_ind_i.at({k}) = + index_f.at({k}) + + dfi::directions[i][k]; + } + + + /* + * ------------------------------------------------ + * Particle momentum projected onto e_i + * ------------------------------------------------ + */ + + auto u_p_e = + saw::math::dot( + e_i, + vel_p_old + ); + + + /* + * ------------------------------------------------ + * Momentum-exchange population + * + * f_i(x) + * + * + f_-i(x + e_i) + * + * with moving-wall correction. + * ------------------------------------------------ + */ + + saw::data dfs_added = + dfs.at({i}) + * ( + saw::data{1.0} + - u_p_e.at({}) + ) + + + dfs_old_f + .at(n_ind_i) + .at({i_opp}) + * ( + saw::data{1.0} + + u_p_e.at({}) + ); + + + /* + * Convert scalar into vector multiplication. + */ + + saw::data> dfs_added_v; + + dfs_added_v.at({}) = + dfs_added; + + + /* + * ------------------------------------------------ + * Delta p_i + * + * ------------------------------------------------ + */ + + momentum = + momentum + + e_i * dfs_added_v; + } + + + /* + * ------------------------------------------------ + * Since both i and i_opp were included above, + * compensate for double counting. + * + * This is the important change from your original + * implementation. + * ------------------------------------------------ + */ + + + + momentum = + momentum * half; + + + /* + * ------------------------------------------------ + * PSM solid fraction + * + * por = 1 -> fluid + * por = 0 -> particle + * + * Therefore: + * + * flip_por = 1 - por + * ------------------------------------------------ + */ + + auto flip_por = + one - por; + + + /* + * ------------------------------------------------ + * Force transferred TO FLUID + * ------------------------------------------------ + */ + + auto& force = + force_f.at(index_f); + + force = + momentum * flip_por; + + + /* + * ------------------------------------------------ + * Force transferred TO PARTICLE + * + * Newton's third law. + * ------------------------------------------------ + */ + + force_p = + force_p - force; + + }, + start, + stop + ); + + + /* + * ------------------------------------------------------------ + * Particle acceleration + * ------------------------------------------------------------ + */ + + auto& pirb_acc = + pirb.template get<"acceleration">(); + + pirb_acc = + force_p + / pg.template get<"total_mass">().at({}); + + + /* + * ------------------------------------------------------------ + * Integrate particle + * ------------------------------------------------------------ + */ + + for(saw::data i{0u}; + i < sub_steps; + ++i){ + + verlet_step_lambda( + pi, + ts + ); + } + } +} }; } diff --git a/modules/core/c++/write_csv.hpp b/modules/core/c++/write_csv.hpp index a60e208..aa28613 100644 --- a/modules/core/c++/write_csv.hpp +++ b/modules/core/c++/write_csv.hpp @@ -91,7 +91,7 @@ struct lbm_csv_writer> { template struct lbm_csv_writer> { static saw::error_or apply(std::ostream& csv_file, const saw::data>& field){ - static_assert(D > 0, "Non-dimensionality is bad for velocity."); + static_assert(D > 0, "Non-dimensionality is bad for momentum."); // csv_file<<"VECTORS "<> { return saw::make_void(); } static saw::error_or apply(std::ostream& vtk_file, const saw::data>& field){ - static_assert(D > 0, "Non-dimensionality is bad for velocity."); + static_assert(D > 0, "Non-dimensionality is bad for momentum."); static_assert(D <= 3, "4th dimension as well. Mostly due to vtk."); // vtk_file<<"VECTORS "<