From 38783c5a86033d63fcbd3af57bb037ff665bfa01 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Tue, 28 Jul 2026 13:48:10 +0200 Subject: Adding fplbm --- default.nix | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/default.nix b/default.nix index f5898ed..241d84a 100644 --- a/default.nix +++ b/default.nix @@ -126,6 +126,11 @@ in rec { inherit pname version stdenv forstio adaptive-cpp; inherit kel; }; + + moving_poiseulle_particles_2d_fplbm_gpu = pkgs.callPackage ./examples/moving_poiseulle_particles_2d_fplbm_gpu/.nix/derivation.nix { + inherit pname version stdenv forstio adaptive-cpp; + inherit kel; + }; poiseulle_particles_2d_psm_gpu = pkgs.callPackage ./examples/poiseulle_particles_2d_psm_gpu/.nix/derivation.nix { inherit pname version stdenv forstio adaptive-cpp; @@ -163,11 +168,6 @@ in rec { ]; }; - poiseulle_moving_particle_2d_psm_gpu = pkgs.callPackage ./examples/poiseulle_moving_particle_2d_psm_gpu/.nix/derivation.nix { - inherit pname version stdenv forstio adaptive-cpp; - inherit kel; - }; - poiseulle_particles_2d_gpu = pkgs.callPackage ./examples/poiseulle_particles_2d_gpu/.nix/derivation.nix { inherit pname version stdenv forstio adaptive-cpp; inherit kel; -- cgit v1.2.3 From 6fa5e2f97c215a0c56679514af74e88f81200d47 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Tue, 28 Jul 2026 20:28:48 +0200 Subject: Fixing up fplbm --- modules/core/c++/fplbm.hpp | 76 ++++++++-------------------------------------- 1 file changed, 13 insertions(+), 63 deletions(-) diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp index 34806aa..31b66c4 100644 --- a/modules/core/c++/fplbm.hpp +++ b/modules/core/c++/fplbm.hpp @@ -114,8 +114,8 @@ private: public: component() = default; - template - void apply(const saw::data& field, const saw::data& macros, saw::data> index, saw::data time_step) const { + 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); @@ -129,76 +129,26 @@ public: // Temporary find a better way auto& force_f = macros.template get<"force">(); - /** * The other methods all work with a flipped porosity. * For compat reasons this is also flipped */ - auto porosity = por_f.at(index); - saw::data> one; - one.at({}) = 1.0; - auto flip_porosity = one - porosity; - - saw::data>& rho = rho_f.at(index); - - saw::data> half; - half.at({}).set(0.5); - saw::data>& vel = vel_f.at(index);// + total_force * ( half / rho ); - - compute_rho_u(dfs_old_f.at(index),rho,vel); - auto eq = equilibrium(rho,vel); - - using dfi = df_info; - - saw::data> min_two; - min_two.at({}).set(-2); - - // Maybe ? - auto& force = force_f.at(index); - force = vel * rho * min_two * flip_porosity; - - saw::data> dfi_inv_cs2; - dfi_inv_cs2.at({}).set(dfi::inv_cs2); - - // auto vel = vel_f.at(index); - - for(uint64_t i = 0u; i < Descriptor::Q; ++i){ - // saw::data ci_min_u{0}; - saw::data> ci; - for(uint64_t d = 0u; d < Descriptor::D; ++d){ - ci.at({{d}}).set(static_cast::type>(dfi::directions[i][d])); - } - auto ci_dot_u = saw::math::dot(ci,vel); - - // saw::data> F_i; - // F_i = f * ((c_i - u) * ics2 + * c_i * ics2 * ics2) * w_i; - saw::data> w; - w.at({}).set(dfi::weights[i]); - - /* - saw::data> F_i_sum; - for(uint64_t d = 0u; d < Descriptor::D; ++d){ - saw::data> F_i_d; - F_i_d.at({}) = F_i.at({{d}}); - F_i_sum = F_i_sum + F_i_d; - } - */ - auto term1 = (ci-vel) * dfi_inv_cs2; - auto term2 = ci * (ci_dot_u * dfi_inv_cs2 * dfi_inv_cs2); - - auto force_projection = saw::math::dot(term1 + term2, force); + // TODO - Change to tuple later + auto parts = pg.template get<"particles">(); - auto F_i = w * force_projection; + auto parts_size = parts.meta.at({0u}); + auto pi = parts.at(index); - dfs.at({i}) = dfs.at({i}) - + frequency_ - * (eq.at(i) - dfs.at({i}) ) - + F_i.at({}) - * (saw::data{1} - saw::data{0.5f} * frequency_); + auto eo_aabb = particle_aabb::calculate(pg,{{0u}},vel_f.meta()); + if(eo_aabb.is_error()){ + return; } - } + auto& aabb = eo_aabb.get_value(); + iterator::apply([&](const auto& index_f) -> void { + }, start, stop); + } }; template -- cgit v1.2.3 From 87b0148bb6814171983b537a456f7b5004833529 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Thu, 30 Jul 2026 20:27:29 +0200 Subject: Ok, probably an init error again --- .../sim.cpp | 65 +++++-------- modules/core/c++/fplbm.hpp | 107 +++++++++++++++------ 2 files changed, 100 insertions(+), 72 deletions(-) diff --git a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp index 3164945..42495c6 100644 --- a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp +++ b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp @@ -48,7 +48,8 @@ template using MacroStruct = Struct< Member, "velocity">, Member, "density">, - Member, "porosity"> + Member, "porosity">, + Member, "force"> >; //template @@ -61,13 +62,12 @@ using ParticleSpheroidGroup = ParticleGroup>; } template -saw::error_or setup_initial_conditions( +saw::error_or init( saw::data>& fields, saw::data>& macros, saw::data>& particles ){ auto& info_f = fields.template get<"info">(); - auto& porous_f = macros.template get<"porosity">(); // Set everything as walls iterator::apply( [&](auto& index){ @@ -107,42 +107,34 @@ saw::error_or setup_initial_conditions( {{0u,1u}} ); // - auto& df_f = fields.template get<"dfs_old">(); + auto& dfs_old_f = fields.template get<"dfs_old">(); + auto& dfs_f = fields.template get<"dfs">(); 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">(); iterator::apply( [&](auto& index){ - auto& df = df_f.at(index); + auto& dfs = dfs_f.at(index); auto& rho = rho_f.at(index); por_f.at(index).at({}) = {1}; rho.at({}) = {1}; auto& vel = vel_f.at(index); - auto eq = equilibrium(rho,vel); - - df = eq; - }, - {},// 0-index - df_f.get_dims() - ); - - iterator::apply( - [&](auto& index){ - auto& df = df_f.at(index); - auto& rho = rho_f.at(index); - rho.at({}) = {1}; - auto& vel = vel_f.at(index); - if(info_f.at(index).get() == 2u){ - vel.at({{0u}}) = 0.0; + for(uint64_t i{0u}; i < Desc::D; ++i){ + vel.at({{i}}).set(0); } auto eq = equilibrium(rho,vel); - df = eq; + dfs = eq; + dfs_old_f.at(index) = dfs; + auto& force = force_f.at(index); + for(uint64_t i{0u}; i < Desc::D; ++i){ + force.at({{i}}).set(0); + } }, {},// 0-index - df_f.get_dims(), - {{1u,1u}} + df_f.get_dims() ); saw::data> rad; @@ -181,19 +173,19 @@ saw::error_or step( auto& porous_f = macros.template get<"porosity">(); q.submit([&](acpp::sycl::handler& h){ - component> hlbm_reset; + component> fplbm_reset; h.parallel_for(acpp::sycl::range{dim_x,dim_y}, [=](acpp::sycl::id idx){ saw::data> index; for(uint64_t i = 0u; i < Desc::D; ++i){ index.at({{i}}).set(idx[i]); } - hlbm_reset.apply(fields,macros,index,t_i); + fplbm_reset.apply(fields,macros,index,t_i); }); }).wait(); q.submit([&](acpp::sycl::handler& h){ - component> hlbm_one_part; + component> fplbm_one_part; h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){ saw::data> index; @@ -201,13 +193,13 @@ saw::error_or step( index.at({{i}}).set(idx[i]); } - hlbm_one_part.apply(fields,macros,particles,index,t_i,{16u}); + fplbm_one_part.apply(fields,macros,particles,index,t_i,{1u}); }); }).wait(); // auto coll_ev = q.submit([&](acpp::sycl::handler& h){ - component> collision{0.8}; + component> collision{0.8}; component> bb; component,encode::Sycl> flow_in{ @@ -221,7 +213,6 @@ saw::error_or step( }; component,encode::Sycl> flow_out{1.0}; - h.parallel_for(acpp::sycl::range{dim_x,dim_y}, [=](acpp::sycl::id idx){ saw::data> index; for(uint64_t i = 0u; i < Desc::D; ++i){ @@ -253,14 +244,6 @@ saw::error_or step( }); }).wait(); - - // Step - /* - q.submit([&](acpp::sycl::handler& h){ - // h.depends_on(collision_ev); - }).wait(); - */ - return saw::make_void(); } } @@ -278,7 +261,7 @@ saw::error_or lbm_main(int argc, char** argv){ } auto& lbm_dir = eo_lbm_dir.get_value(); - auto out_dir = lbm_dir / "moving_poiseulle_particles_2d_hlbm_gpu"; + auto out_dir = lbm_dir / "moving_poiseulle_particles_2d_fplbm_gpu"; { std::error_code ec; @@ -327,7 +310,7 @@ saw::error_or lbm_main(int argc, char** argv){ sycl_q.wait(); { - auto eov = setup_initial_conditions(*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_parts_data_ptr); + auto eov = init(*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_parts_data_ptr); if(eov.is_error()){ return eov; } @@ -421,6 +404,8 @@ saw::error_or lbm_main(int argc, char** argv){ } }); }).wait(); + + // EPoll wait.poll(); if(print_status){ std::cout<<"Status: "<().get() * 100 / time_steps.get())<<"%"<(); + auto& force = force_f.at(index); + for(uint64_t i{0u}; i < Descriptor::D; ++i){ + force.at({{i}}) = 0.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& 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); } }; @@ -52,19 +73,18 @@ public: template void apply(const saw::data& field, const saw::data& macros, saw::data> index, saw::data time_step) 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">(); - saw::data>& rho = rho_f.at(index); + + auto& vel_f = macros.template get<"velocity">(); saw::data>& vel = vel_f.at(index); - compute_rho_u(dfs_old_f.at(index),rho,vel); + compute_rho_u(dfs,rho,vel); auto eq = equilibrium(rho,vel); using dfi = df_info; @@ -72,36 +92,34 @@ public: auto& force_f = macros.template get<"force">(); auto& force = force_f.at(index); - auto& porosity_f = macros.template get<"porosity">(); - auto& porosity = porosity_f.at(index); + auto& por_f = macros.template get<"porosity">(); + auto& por = por_f.at(index); saw::data> dfi_inv_cs2; dfi_inv_cs2.at({}).set(dfi::inv_cs2); - for(uint64_t i = 0u; i < Descriptor::Q; ++i){ - // saw::data ci_min_u{0}; + for(uint64_t i{0u}; i < Descriptor::Q; ++i){ saw::data> ci; - for(uint64_t d = 0u; d < Descriptor::D; ++d){ + for(uint64_t d{0u}; d < Descriptor::D; ++d){ ci.at({{d}}).set(static_cast::type>(dfi::directions[i][d])); } auto ci_dot_u = saw::math::dot(ci,vel); - // saw::data> F_i; - // F_i = f * (c_i - u * ics2 + * c_i * ics2 * ics2) * w_i; saw::data> w; w.at({}).set(dfi::weights[i]); - auto F_i_d = saw::math::dot(force * w, (ci - vel * dfi_inv_cs2 + ci * ci_dot_u * dfi_inv_cs2 * dfi_inv_cs2 )); - /* - saw::data> F_i_sum; - for(uint64_t d = 0u; d < Descriptor::D; ++d){ - saw::data> F_i_d; - F_i_d.at({}) = F_i.at({{d}}); - F_i_sum = F_i_sum + F_i_d; - } - */ + auto term1 = (ci-vel) * dfi_inv_cs2; + auto term2 = ci * (ci_dot_u * dfi_inv_cs2 * dfi_inv_cs2); + + auto force_projection = saw::math::dot(term1 + term2, force); + + auto F_i = w * force_projection; - dfs.at({i}) = dfs.at({i}) + frequency_ * (eq.at(i) - dfs.at({i}) ) + F_i_d.at({}) * (saw::data{1} - saw::data{0.5f} * frequency_); + dfs.at({i}) = dfs.at({i}) + + frequency_ + * (eq.at(i) - dfs.at({i}) ) + + F_i.at({}) + * (saw::data{1} - saw::data{0.5f} * frequency_); } } }; @@ -119,14 +137,12 @@ public: // 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& 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">(); - - // Temporary find a better way auto& force_f = macros.template get<"force">(); /** @@ -136,8 +152,19 @@ public: // TODO - Change to tuple later auto parts = pg.template get<"particles">(); - auto parts_size = parts.meta.at({0u}); - auto pi = parts.at(index); + 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">(); + + saw::data> sss; + sss.at({}) = sub_steps.template cast_to(); + auto vel_s = (pirb_pos-pirb_pos_old) / sss; + + 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()){ @@ -145,9 +172,28 @@ public: } auto& aabb = eo_aabb.get_value(); - iterator::apply([&](const auto& index_f) -> void { + saw::data> two; + two.at({}).set(2); - }, start, stop); + saw::data> one; + one.at({}) = 1.0; + + saw::data> eps; + eps.at({}) = 1.5f; + + 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; + + force = ( vel_s * rho - vel * rho ) * two * flip_por; + }, aabb.template get<"a">(), aabb.template get<"b">()); } }; @@ -208,15 +254,12 @@ public: saw::data> min_two; min_two.at({}).set(-2); - // Maybe ? auto& force = force_f.at(index); force = vel * rho * min_two * flip_porosity; saw::data> dfi_inv_cs2; dfi_inv_cs2.at({}).set(dfi::inv_cs2); - // auto vel = vel_f.at(index); - for(uint64_t i = 0u; i < Descriptor::Q; ++i){ // saw::data ci_min_u{0}; saw::data> ci; -- cgit v1.2.3 From ea3564475a0e8431c553fb34d4b680e05c451c8a Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Thu, 30 Jul 2026 21:27:06 +0200 Subject: Turns out, you need to reference the data type you retrieve ._. --- examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp | 4 ++-- examples/poiseulle_particles_2d_gpu/sim.cpp | 4 +++- modules/core/c++/fplbm.hpp | 16 +++++++++++++--- 3 files changed, 18 insertions(+), 6 deletions(-) diff --git a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp index 42495c6..f6c09e9 100644 --- a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp +++ b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp @@ -134,7 +134,7 @@ saw::error_or init( } }, {},// 0-index - df_f.get_dims() + dfs_f.get_dims() ); saw::data> rad; @@ -350,7 +350,7 @@ saw::error_or lbm_main(int argc, char** argv){ auto lsdm_view = make_view(lbm_sycl_macro_data); auto lsdp_view = make_view(lbm_sycl_parts_data); - saw::data time_steps{16u*4096ul}; + saw::data time_steps{32u*2048ul}; auto& info_f = lsd_view.template get<"info">(); for(saw::data i{0u}; i < time_steps and krun; ++i){ diff --git a/examples/poiseulle_particles_2d_gpu/sim.cpp b/examples/poiseulle_particles_2d_gpu/sim.cpp index fd6cdca..f5b49e9 100644 --- a/examples/poiseulle_particles_2d_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_gpu/sim.cpp @@ -25,6 +25,8 @@ template saw::error_or lbm_main(const saw::data& args){ using namespace kel::lbm; + auto& an = args.template get<"args">(); + using dfi = df_info; auto eo_lbm_dir = output_directory(); @@ -33,7 +35,7 @@ saw::error_or lbm_main(const saw::data& args){ } auto& lbm_dir = eo_lbm_dir.get_value(); - auto out_dir = lbm_dir / "poiseulle_particles_2d_gpu" / "hlbm"; + auto out_dir = lbm_dir / "poiseulle_particles_2d_gpu" / an.template get<"coupling">().stl_view(); { std::error_code ec; diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp index 3c8e173..cfb0062 100644 --- a/modules/core/c++/fplbm.hpp +++ b/modules/core/c++/fplbm.hpp @@ -1,8 +1,7 @@ #pragma once #include "common.hpp" -#include "particle/aabb.hpp" -#include "particle/porosity.hpp" +#include "particle/particle.hpp" #include "iterator.hpp" namespace kel { @@ -153,7 +152,7 @@ public: auto parts = pg.template get<"particles">(); auto parts_size = parts.meta().at({0u}); - auto pi = parts.at(index); + auto& pi = parts.at(index); auto& pirb = pi.template get<"rigid_body">(); auto& pirb_pos = pirb.template get<"position">(); @@ -181,6 +180,8 @@ public: saw::data> eps; eps.at({}) = 1.5f; + 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); @@ -193,7 +194,16 @@ public: auto flip_por = one - por; 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); + } } }; -- cgit v1.2.3