From c0549d71b2109f10c1238db8b22362e7826ba61b Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Sun, 5 Jul 2026 15:59:23 +0200 Subject: Just rename from lib to modules --- modules/core/c++/hlbm.hpp | 168 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 168 insertions(+) create mode 100644 modules/core/c++/hlbm.hpp (limited to 'modules/core/c++/hlbm.hpp') diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp new file mode 100644 index 0000000..18c32a5 --- /dev/null +++ b/modules/core/c++/hlbm.hpp @@ -0,0 +1,168 @@ +#pragma once + +#include "macroscopic.hpp" +#include "component.hpp" +#include "equilibrium.hpp" + +#include "particle/particle.hpp" + +#include + +namespace kel { +namespace lbm { +namespace cmpt { +struct HlbmInit {}; +struct Hlbm {}; +struct HlbmParticle {}; + +struct HlbmOneParticleMomentumExchange {}; +} + +template +class component final { +private: + typename saw::native_data_type::type relaxation_; + saw::data frequency_; +public: + component(typename saw::native_data_type::type relaxation__): + relaxation_{relaxation__}, + frequency_{typename saw::native_data_type::type(1) / relaxation_} + {} + + 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& particle_N_f = field.template get<"particle_N">(); + auto& particle_D_f = field.template get<"particle_D">(); + + auto& por = porosity_f.at(index); + por = {}; + + auto& pnf = particle_N_f.at(index); + pnf = {}; + + auto& pnd = particle_D_f.at(index); + pnd = {}; + } +}; + +/** + * HLBM collision operator for LBM + */ +template +class component final { +private: + typename saw::native_data_type::type relaxation_; + saw::data frequency_; +public: + component(typename saw::native_data_type::type relaxation__): + relaxation_{relaxation__}, + frequency_{typename saw::native_data_type::type(1) / relaxation_} + {} + + template + void apply(const saw::data& field, const saw::data& macros, saw::data> index, saw::data time_step) const { + + bool is_even = ((time_step.get() % 2) == 0); + + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + auto& particle_N_f = field.template get<"particle_N">(); + auto& particle_D_f = field.template get<"particle_D">(); + + auto& porosity_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& porosity = porosity_f.at(index); + saw::data> one; + one.at({}) = 1.0; + auto flip_porosity = one - porosity; + + auto& N = particle_N_f.at(index); + auto& D = particle_D_f.at(index); + // Convex combination of velocities + vel = vel * porosity + [&]() -> saw::data> { + return (D.at({}).get() > 0.0 ? N * flip_porosity / D : N); + }(); + // Equilibrium + auto eq = equilibrium(rho,vel); + + for(uint64_t i = 0u; i < Desc::Q; ++i){ + 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})); + } + + porosity.at({}) = 1.0; + D.at({}) = 0.0; + N = {}; + } +}; + +namespace impl { + +} + +template +class component final { +private: + template + void apply_i(const saw::data& field, const saw::data& macros, const saw::data& part_groups, saw::data> index, saw::data time_step) const { + + } +public: + /* + template + void apply_i() + */ + + template + void apply(const saw::data& field, const saw::data& macros, const saw::data& part_group, saw::data> index, saw::data time_step) const { + /// Figure out how to access the particle list + // auto& p = particles.at(i); + + /// Iterate over the grid bounds + // auto& grid = p.template get<"grid">(); + + auto& part_spheroid_group = part_group; + auto& mvel = macros.template get<"velocity">(); + { + auto& parts = part_spheroid_group.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">(); + + saw::data> start; + saw::data> stop; + + auto aabb = particle_aabb::calculate(part_spheroid_group,index,mvel.meta()); + /// Ok, I iterate over the space which covers our particle? So lower bounds to upper bounds + + iterator::apply([&](const auto& index){ + // ask for the d_k value here. + // For every value im iterating over I need sth + std::cout<<"Pos: "< +class component final { +public: + template + void apply(const saw::data& field, const saw::data& macros, const saw::data& part_group, saw::data> index, saw::data time_step) const { + // + } +}; +} +} -- cgit v1.2.3