From eb97ceef955fa1c5e3794c74fc9894fda1ce6f21 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Tue, 2 Jun 2026 18:15:13 +0200 Subject: Doing aabb computations --- default.nix | 14 + .../.nix/derivation.nix | 41 ++ examples/settling_spheres_2d_hlbm_gpu/SConscript | 34 ++ examples/settling_spheres_2d_hlbm_gpu/SConstruct | 81 ++++ examples/settling_spheres_2d_hlbm_gpu/sim.cpp | 469 +++++++++++++++++++++ lib/core/c++/chunk.hpp | 5 + lib/core/c++/hlbm.hpp | 43 +- lib/core/c++/particle/aabb.hpp | 49 +++ lib/core/c++/particle/porosity.hpp | 57 ++- 9 files changed, 778 insertions(+), 15 deletions(-) create mode 100644 examples/settling_spheres_2d_hlbm_gpu/.nix/derivation.nix create mode 100644 examples/settling_spheres_2d_hlbm_gpu/SConscript create mode 100644 examples/settling_spheres_2d_hlbm_gpu/SConstruct create mode 100644 examples/settling_spheres_2d_hlbm_gpu/sim.cpp create mode 100644 lib/core/c++/particle/aabb.hpp diff --git a/default.nix b/default.nix index d233d90..c8e99f5 100644 --- a/default.nix +++ b/default.nix @@ -161,11 +161,23 @@ in rec { inherit pname version stdenv forstio adaptive-cpp; inherit kel; }; + + settling_cubes_2d_hlbm_gpu = pkgs.callPackage ./examples/settling_spheres_2d_hlbm_gpu/.nix/derivation.nix { + inherit pname version stdenv forstio adaptive-cpp; + inherit kel; + }; heterogeneous_computing = pkgs.callPackage ./examples/heterogeneous_computing/.nix/derivation.nix { inherit pname version stdenv forstio adaptive-cpp; inherit kel; }; + + settling_particles = pkgs.symlinkJoin { + name = "kel-lbm-settling_particles-${version}"; + paths = [ + examples.settling_cubes_2d_ibm_gpu + ]; + }; }; debug = { @@ -190,6 +202,8 @@ in rec { examples.poiseulle_particles_2d_ibm_gpu ]; }; + + }; default = release.examples; diff --git a/examples/settling_spheres_2d_hlbm_gpu/.nix/derivation.nix b/examples/settling_spheres_2d_hlbm_gpu/.nix/derivation.nix new file mode 100644 index 0000000..26b41f2 --- /dev/null +++ b/examples/settling_spheres_2d_hlbm_gpu/.nix/derivation.nix @@ -0,0 +1,41 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, python3 +, pname +, version +, adaptive-cpp +, kel +}: + +stdenv.mkDerivation { + pname = pname + "-examples-" + "setting_spheres_2d_hlbm_gpu"; + inherit version; + src = ./..; + + nativeBuildInputs = [ + scons + clang-tools + python3 + ]; + + buildInputs = [ + forstio.core + forstio.async + forstio.codec + forstio.codec-unit + forstio.io + forstio.remote + forstio.remote-filesystem + forstio.codec-json + adaptive-cpp + kel.lbm.core + kel.lbm.sycl + ]; + + preferLocalBuild = true; + + outputs = [ "out" "dev" ]; +} diff --git a/examples/settling_spheres_2d_hlbm_gpu/SConscript b/examples/settling_spheres_2d_hlbm_gpu/SConscript new file mode 100644 index 0000000..ae08056 --- /dev/null +++ b/examples/settling_spheres_2d_hlbm_gpu/SConscript @@ -0,0 +1,34 @@ +#!/bin/false + +import os +import os.path +import glob + + +Import('env') + +dir_path = Dir('.').abspath + +# Environment for base library +examples_env = env.Clone(); +examples_env['CXX'] = 'syclcc-clang'; +examples_env['CXXFLAGS'] += ['-O3']; + +examples_env.sources = sorted(glob.glob(dir_path + "/*.cpp")) +examples_env.headers = sorted(glob.glob(dir_path + "/*.hpp")) + +env.sources += examples_env.sources; +env.headers += examples_env.headers; + +# Cavity2D +examples_objects = []; +examples_env.add_source_files(examples_objects, ['sim.cpp'], shared=False); +examples_env.poiseulle_2d_gpu = examples_env.Program('#bin/settling_cubes_2d_ibm_gpu', [examples_objects]); + +# Set Alias +env.examples = [ + examples_env.poiseulle_2d_gpu +]; +env.Alias('examples', env.examples); +env.targets += ['examples']; +env.Install('$prefix/bin/', env.examples); diff --git a/examples/settling_spheres_2d_hlbm_gpu/SConstruct b/examples/settling_spheres_2d_hlbm_gpu/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/settling_spheres_2d_hlbm_gpu/SConstruct @@ -0,0 +1,81 @@ +#!/usr/bin/env python3 + +import sys +import os +import os.path +import glob +import re + + +if sys.version_info < (3,): + def isbasestring(s): + return isinstance(s,basestring) +else: + def isbasestring(s): + return isinstance(s, (str,bytes)) + +def add_kel_source_files(self, sources, filetype, lib_env=None, shared=False, target_post=""): + + if isbasestring(filetype): + dir_path = self.Dir('.').abspath + filetype = sorted(glob.glob(dir_path+"/"+filetype)) + + for path in filetype: + target_name = re.sub( r'(.*?)(\.cpp|\.c\+\+)', r'\1' + target_post, path ) + if shared: + target_name+='.os' + sources.append( self.SharedObject( target=target_name, source=path ) ) + else: + target_name+='.o' + sources.append( self.StaticObject( target=target_name, source=path ) ) + pass + +def isAbsolutePath(key, dirname, env): + assert os.path.isabs(dirname), "%r must have absolute path syntax" % (key,) + +env_vars = Variables( + args=ARGUMENTS +) + +env_vars.Add('prefix', + help='Installation target location of build results and headers', + default='/usr/local/', + validator=isAbsolutePath +) + +env_vars.Add('build_examples', + help='If examples should be built', + default="true" +) + +env=Environment(ENV=os.environ, variables=env_vars, CPPPATH=[], + CPPDEFINES=['SAW_UNIX'], + CXXFLAGS=[ + '-std=c++20', + '-g', + '-Wall', + '-Wextra' + ], + LIBS=[ + 'forstio-core', + 'forstio-async', + 'forstio-io' + ] +); +env.__class__.add_source_files = add_kel_source_files +env.Tool('compilation_db'); +env.cdb = env.CompilationDatabase('compile_commands.json'); + +env.objects = []; +env.sources = []; +env.headers = []; +env.targets = []; + +Export('env') +SConscript('SConscript') + +env.Alias('cdb', env.cdb); +env.Alias('all', [env.targets]); +env.Default('all'); + +env.Alias('install', '$prefix') diff --git a/examples/settling_spheres_2d_hlbm_gpu/sim.cpp b/examples/settling_spheres_2d_hlbm_gpu/sim.cpp new file mode 100644 index 0000000..eb3f377 --- /dev/null +++ b/examples/settling_spheres_2d_hlbm_gpu/sim.cpp @@ -0,0 +1,469 @@ +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace kel { +namespace lbm { + +constexpr uint64_t dim_x = 512ul; +constexpr uint64_t dim_y = dim_x * 2ul; + +constexpr uint64_t particle_amount = 16ul; + +namespace sch { +using namespace saw::schema; + +using InfoChunk = Chunk; + +template +using DfChunk = Chunk, 1u, dim_x, dim_y>; + +template +using ScalarChunk = Chunk, 0u, dim_x, dim_y>; + +template +using VectorChunk = Chunk, 0u, dim_x, dim_y>; + +template +using ChunkStruct = Struct< + Member, + Member, "dfs">, + Member, "dfs_old">, + Member, "p_N">, + Member, "p_D"> +>; + +template +using VelChunk = Chunk, 0u, dim_x, dim_y>; + +template +using RhoChunk = Chunk, 0u, dim_x, dim_y>; + +template +using PorChunk = Chunk, 0u, dim_x, dim_y>; + +template +using MacroStruct = Struct< + Member, "velocity">, + Member, "density">, + Member, "porosity"> +>; + +template +using ParticleSpheroidGroup = ParticleGroup< + T, + Desc::D, + sch::ParticleCollisionSpheroid +>; + +template +using ParticleGroups = Tuple< + ParticleSpheroidGroup +>; + + +} + +template +saw::error_or setup_initial_conditions( + saw::data>& fields, + saw::data>& macros, + saw::data>& particles +){ + auto& info_f = fields.template get<"info">(); + // Set everything as walls + iterator::apply( + [&](auto& index){ + info_f.at(index).set(1u); + }, + {}, + info_f.get_dims(), + {} + ); + // Fluid + iterator::apply( + [&](auto& index){ + info_f.at(index).set(2u); + }, + {}, + info_f.get_dims(), + {{1u,1u}} + ); + // + auto& df_f = fields.template gStarted hearing about similar cases not long after, many of which were successful on the part of the crooks. Scary stuff. et<"dfs_old">(); + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"velocity">(); + + 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); + auto eq = equilibrium(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims() + ); + + // Particles + { + saw::data> dense_p; + dense_p.at({}).set(1); + // auto& spheroid_group = particles.template get<0u>(); + auto& spheroid_group = particles; + + spheroid_group = create_spheroid_particle_group( + dense_p, + {64u} + ); + + { + auto& p = spheroid_group.template get<"particles">(); + + /// 16 if I remember correctly? + p = {{{particle_amount}}}; + + iterator<1u>::apply( + [&](auto& index){ + // Set Pos here? + auto& p_ind = p.at(index); + + auto& p_rb = p_ind.template get<"rigid_body">(); + auto& p_pos = p_rb.template get<"position">(); + + saw::data> p; + uint64_t x = index.at({0u}).get() % 8u; + uint64_t y = index.at({1u}).get() / 8u; + p.at({{0u}}).set( 0.2 + 0.6 / (x-1u) ); + p.at({{1u}}).set( 0.8 + y * 0.1 ); + + p_pos = p; + + auto& p_pos_old = p_rb.template get<"position_old">(); + p_pos_old = p_pos; + }, + {}, + p.meta() + ); + } + } + // Particle in hacky flavour + {} + + return saw::make_void(); +} + +template +saw::error_or step( + saw::data>,encode::Sycl>& fields, + saw::data>,encode::Sycl>& macros, + saw::data>,encode::Sycl>& p_group, + saw::data t_i, + device& dev +){ + auto& q = dev.get_handle(); + auto& info_f = fields.template get<"info">(); + + auto& parts = p_group.template get<"particles">(); + auto& p_mask = p_group.template get<"mask">(); + auto& vels = macros.template get<"velocity">(); + auto& forces = macros.template get<"force">(); + + auto p_meta = parts.meta(); + q.submit([&](acpp::sycl::handler& h){ + + 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]); + } + + // Reset the force to zero + forces.at(index) = {}; + }); + }).wait(); + + q.submit([&](acpp::sycl::handler& h){ + h.parallel_for(acpp::sycl::range<1u>{p_meta.at({0u}).get()}, [=](acpp::sycl::id<1u> idx){ + + saw::data> index; + for(uint64_t i = 0u; i < 1u; ++i){ + index.at({{i}}).set(idx[i]); + } + + auto& pi = parts.at(index); + auto& pirb = pi.template get<"rigid_body">(); + + auto& p_pos = pirb.template get<"position">(); + auto& p_rot = pirb.template get<"rotation">(); + auto& p_acc = pirb.template get<"acceleration">(); + + + + saw::data> delta_t; + delta_t.at({}).set(1.0f); + verlet_step_lambda(p,delta_t); + }); + }).wait(); + + // auto coll_ev = + q.submit([&](acpp::sycl::handler& h){ + // Need nicer things to handle the flow. I see improvement here + component> collision{0.8}; + component> bb; + + 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]); + } + + auto info = info_f.at(index); + + switch(info.get()){ + case 0u: + break; + case 1u: + bb.apply(fields,index,t_i); + break; + case 2u: + // collision.apply(fields,macros,index,t_i); + break; + default: + break; + } + }); + }).wait(); + + // Step + /* + q.submit([&](acpp::sycl::handler& h){ + // h.depends_on(collision_ev); + }).wait(); + */ + + return saw::make_void(); +} +} +} + +template +saw::error_or lbm_main(int argc, char** argv){ + using namespace kel::lbm; + + using dfi = df_info; + + auto eo_lbm_env = lbm_directory(); + if(eo_lbm_env.is_error()){ + return std::move(eo_lbm_env.get_error()); + } + auto& lbm_env = eo_lbm_env.get_value(); + + auto out_dir = lbm_env.data_dir / "settling_cubes_2d_ibm_gpu"; + + { + std::error_code ec; + std::filesystem::create_directories(out_dir,ec); + if(ec != std::errc{}){ + return saw::make_error("Could not create output directory"); + } + } + + converter conv { + // delta_x + {{1.0}}, + // delta_t + {{1.0}} + }; + + print_lbm_meta(conv,{0.05},{0.01},{0.4 * dim_y}); + + // saw::data> meta{{dim_x,dim_y}}; + auto lbm_data_ptr = saw::heap>>(); + auto lbm_macro_data_ptr = saw::heap>>(); + auto lbm_particle_data_ptr = saw::heap>>(); + // auto lbm_particles_ptr = saw::heap,part_count>>>(); + // saw::data> p_mask; + + std::cout<<"Estimated Bytes of LBM Fields: "<,sch::MacroStruct>().get()<(*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_vtk_file(out_dir,"initial_state",0u,*lbm_data_ptr); + if(eov.is_error()){ + return eov; + } + } + + saw::data, encode::Sycl> lbm_sycl_data{sycl_q}; + saw::data, encode::Sycl> lbm_sycl_macro_data{sycl_q}; + saw::data, encode::Sycl> lbm_sycl_particle_data{sycl_q}; + + { + auto eov = dev.copy_to_device(*lbm_data_ptr,lbm_sycl_data); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = dev.copy_to_device(*lbm_macro_data_ptr,lbm_sycl_macro_data); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = dev.copy_to_device(*lbm_particle_data_ptr,lbm_sycl_particle_data); + if(eov.is_error()){ + return eov; + } + } + sycl_q.wait(); + + auto lsd_view = make_view(lbm_sycl_data); + auto lsdm_view = make_view(lbm_sycl_macro_data); + auto lsdp_view = make_view(lbm_sycl_particle_data); + + { + auto eov = write_vtk_file(out_dir,"ms",0u, *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + saw::data time_steps{16u*4096ul}; + + auto& info_f = lsd_view.template get<"info">(); + + for(saw::data i{0u}; i < time_steps and krun; ++i){ + // BC + Collision + { + auto eov = step(lsd_view,lsdm_view,lsdp_view,i,dev); + if(eov.is_error()){ + return eov; + } + } + sycl_q.wait(); + if(i.get()%1u ==0u){ + { + auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_vtk_file(out_dir,"m",i.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + } + /* + { + { + auto eov = dev.copy_to_host(lbm_sycl_data,*lbm_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_csv_file(out_dir,"lbm",i.get(), *lbm_data_ptr); + if(eov.is_error()){ + return eov; + } + } + } + */ + // Stream + sycl_q.submit([&](acpp::sycl::handler& h){ + component> stream; + + 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]); + } + + auto info = info_f.at(index); + + if(info.get() > 0u){ + stream.apply(lsd_view,index,i); + } + }); + }).wait(); + wait.poll(); + if(print_status){ + std::cout<<"Status: "<().get() * 100 / time_steps.get())<<"%"<(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + std::cerr<<"[Error] "< 0u){ + std::cerr<<" - "<& neighbour_at(const data>& index) const { + } + */ + static constexpr auto get_dims(){ return data,Encode>{{Sides...}}; } diff --git a/lib/core/c++/hlbm.hpp b/lib/core/c++/hlbm.hpp index 7590cc2..ed49353 100644 --- a/lib/core/c++/hlbm.hpp +++ b/lib/core/c++/hlbm.hpp @@ -96,20 +96,55 @@ public: } }; -template +namespace impl { + +} + +template class component final { private: - typename saw::native_data_type::type relaxation_; - saw::data frequency_; + 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& particles, saw::data> index, saw::data time_step) const { + void apply(const saw::data& field, const saw::data& macros, const saw::data& part_groups, 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_groups.template get<0>(); + { + auto& parts = part_spheroid_group.template get<"particles">(); + auto parts_size = parts.size(); + + 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; + + /// Ok, I iterate over the space which covers our particle? So lower bounds to upper bounds + for(uint64_t i{0u}; i < Desc::D; ++i){ + + } + + iterator::apply([&](const auto& index){ + // ask for the d_k value here. + // For every value im iterating over I need sth + },start,stop); + + // Check + } } diff --git a/lib/core/c++/particle/aabb.hpp b/lib/core/c++/particle/aabb.hpp new file mode 100644 index 0000000..aec95ca --- /dev/null +++ b/lib/core/c++/particle/aabb.hpp @@ -0,0 +1,49 @@ +#pragma once + +#include "particle.hpp" + +namespace kel { +namespace lbm { +template +class particle_aabb final { +}; + +template::type radius> +class particle_aabb > > final { +public: + using Schema = sch::ParticleGroup>; + + using AABB = Struct< + Member"a">, + Member"b"> + >; +public: + static constexpr saw::data get(const saw::data& p_grp, const saw::data>& i, const saw::data>& meta){ + auto& parts = p_grp.template get<"particles">(); + + auto& pi = parts.at(i); + auto& pirb = pi.template get<"rigid_body">(); + auto& pirb_pos = pirb.template get<"position">(); + + saw::data aabb; + auto& a = aabb.template get<"a">(); + auto& b = aabb.template get<"b">(); + + saw::data> rad_d; + rad_d.at({}).set(radius); + + saw::data> lower; + saw::data> upper; + + for(uint64_t i{0u}; i < D; ++i){ + lower.at({{i}}) = pirb_pos.at({{i}}) >= rad_d.at({}) ? (pirb_pos.at({{i}}) - rad_d.at({})) : saw::data{0}; + a.at({i}) = lower.at({{i}}).template cast_to(); + upper.at({{i}}) = pirb_pos.at({{i}}) + rad_d.at({}); + b.at({i}) = (upper.at({{i}})+saw::data{1}).template cast_to() + } + + return aabb; + } +}; +} +} diff --git a/lib/core/c++/particle/porosity.hpp b/lib/core/c++/particle/porosity.hpp index aa1ce5b..39d9652 100644 --- a/lib/core/c++/particle/porosity.hpp +++ b/lib/core/c++/particle/porosity.hpp @@ -8,7 +8,9 @@ namespace lbm { template class particle_porosity { public: - saw::data> calculate(const saw::data<>& part_group, uint64_t p_i, const saw::data>& lbm_pos){ + + + saw::data> calculate(const saw::data>& part_group, uint64_t p_i, const saw::data>& lbm_pos){ auto& mask = part_group.template get<"mask">(); auto& particles = part_group.template get<"particles">(); @@ -17,7 +19,7 @@ public: auto& part_i_rb = part_i.template get<"rigid_body">(); auto& pirb = part_i_rb.template get<"position">(); - auto& dist = lbm_pos = lbm_pos - pirb; + auto dist = lbm_pos - pirb; // index 0 is at @@ -26,25 +28,58 @@ public: }; -template::type radius> -class particle_porosity> final { +template::type radius, typename saw::native_data_type::type eps> +class particle_porosity> final { public: - saw::data> calculate(const saw::data&, uint64_t i, const saw::data>& lbm_pos){ + saw::data> calculate(const saw::data > >& part_group, uint64_t i, const saw::data>& lbm_pos){ saw::data> por; por.at({}); - saw::data> dps_2; + auto& parts = part_group.template get<"particles">(); + auto& pi = parts.at({i}); + auto& pi_rb = pi.template get<"rigid_body">(); + + auto& pi_rb_pos = pi_rb.template get<"position">(); + + // Basically the queried position minus the center of the particle + auto dist = lbm_pos - pi_rb_pos; + + saw::data> dist_len; for(uint64_t i{0u}; i < D; ++i){ - auto& dps_i = lbm_pos.at({{i}}); - dps_2.at({}) = dps_i * dps_i; + auto& d_i = dist.at({{i}}); + dist_len.at({}) = d_i * d_i; } + dist_len.at({}).set(std::sqrt(dist_len.at({}).get()); - saw::data> rad_2; - rad_2.at({}).set(radius*radius); + saw::data> rad_d{radius}; + saw::data> eps_half{eps *0.5}; + + saw::data> rad_d_eps_p = rad_d + eps_half; + saw::data> rad_d_eps_n = rad_d - eps_half; + + // saw::data> rad_d_eps_p_2 = rad_d_eps_p * rad_d_eps_p; + // saw::data> rad_d_eps_n_2 = rad_d_eps_n * rad_d_eps_n; + // saw::data> rad_2; + // rad_2 = rad_d * rad_d; saw::data> inside; - if(dps_2.at({}).get() < rad_2.at({}).get()){ + if(dist_len.at({}) <= rad_d_eps_n.at({})){ inside.at({}).set(1); + return inside; + } + if(dist_len.at({}) > rad_d_eps_p.at({})){ + inside.at({}).set(0); + return inside + } + + /* + * cos^2 ( ||x-X(t)||_2 - (R-eps/2) ) + */ + { + typename saw::native_data_type::type inner = (std::pi / ( 2.0 * eps )) * (dist_len - rad_d_eps_n).at({}).get(); + auto cos_inner = std::cos(inner); + auto cos_i_2 = cos_inner * cos_inner; + inside.at({}).set(cos_i_2); } return inside; } -- cgit v1.2.3 From d7c88d561d4e2b7a2706d14c11ea5fa4cd6ceaf3 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Tue, 2 Jun 2026 19:16:37 +0200 Subject: Fixing Desc name --- lib/core/c++/hlbm.hpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/lib/core/c++/hlbm.hpp b/lib/core/c++/hlbm.hpp index ed49353..6ae7d80 100644 --- a/lib/core/c++/hlbm.hpp +++ b/lib/core/c++/hlbm.hpp @@ -43,8 +43,8 @@ public: /** * HLBM collision operator for LBM */ -template -class component final { +template +class component final { private: typename saw::native_data_type::type relaxation_; saw::data frequency_; @@ -55,7 +55,7 @@ public: {} template - void apply(const saw::data& field, const saw::data& macros, saw::data> index, saw::data time_step) const { + 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); @@ -68,9 +68,9 @@ public: auto& vel_f = macros.template get<"velocity">(); saw::data>& rho = rho_f.at(index); - saw::data>& vel = vel_f.at(index); + saw::data>& vel = vel_f.at(index); - compute_rho_u(dfs_old_f.at(index), rho, vel); + compute_rho_u(dfs_old_f.at(index), rho, vel); auto& porosity = porosity_f.at(index); saw::data> one; @@ -80,13 +80,13 @@ public: auto& N = particle_N_f.at(index); auto& D = particle_D_f.at(index); // Convex combination of velocities - vel = vel * porosity + [&]() -> saw::data> { + vel = vel * porosity + [&]() -> saw::data> { return (D.at({}).get() > 0.0 ? N * flip_porosity / D : N); }(); // Equilibrium - auto eq = equilibrium(rho,vel); + auto eq = equilibrium(rho,vel); - for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + 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})); } @@ -101,7 +101,7 @@ namespace impl { } template -class component final { +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 { -- cgit v1.2.3 From 1c9d05a7cf13b4bd4ebe5583d10cae85d507ea9c Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Wed, 3 Jun 2026 09:02:33 +0200 Subject: Adding reset of particle after collision. --- examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp index 9375078..b09aa64 100644 --- a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp @@ -13,7 +13,7 @@ namespace lbm { constexpr uint64_t dim_y = 256ul; constexpr uint64_t dim_x = dim_y * 20ul; -constexpr uint64_t particle_amount = 1ul; +// constexpr uint64_t particle_amount = 1ul; namespace sch { using namespace saw::schema; @@ -215,6 +215,21 @@ saw::error_or step( break; case 2u: collision.apply(fields,macros,index,t_i); + { + saw::data> middle; + middle.at({{0u}}) = dim_x * 0.25; + middle.at({{1u}}) = dim_y * 0.5; + + saw::data> ind_vec; + ind_vec.at({{0u}}) = index.at({{0u}}).template cast_to(); + ind_vec.at({{1u}}) = index.at({{1u}}).template cast_to(); + + auto dist = middle - ind_vec; + auto dist_2 = saw::math::dot(dist,dist); + if(dist_2.at({}).get() < dim_y*dim_y*0.01){ + porous_f.at(index).at({}) = 0.0; + } + } break; case 3u: flow_in.apply(fields,index,t_i); -- cgit v1.2.3 From 3700a2bdc86a5c56a4051bc5262878e0858306e8 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Wed, 3 Jun 2026 20:23:27 +0200 Subject: Dangling --- default.nix | 7 +- .../.nix/derivation.nix | 41 ++ .../SConscript | 34 ++ .../SConstruct | 81 ++++ .../moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp | 428 +++++++++++++++++++++ examples/poiseulle_particles_2d_ibm_gpu/sim.cpp | 2 +- examples/settling_spheres_2d_hlbm_gpu/sim.cpp | 10 +- 7 files changed, 594 insertions(+), 9 deletions(-) create mode 100644 examples/moving_poiseulle_particles_2d_hlbm_gpu/.nix/derivation.nix create mode 100644 examples/moving_poiseulle_particles_2d_hlbm_gpu/SConscript create mode 100644 examples/moving_poiseulle_particles_2d_hlbm_gpu/SConstruct create mode 100644 examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp diff --git a/default.nix b/default.nix index c8e99f5..99e306a 100644 --- a/default.nix +++ b/default.nix @@ -122,6 +122,11 @@ in rec { inherit kel; }; + moving_poiseulle_particles_2d_hlbm_gpu = pkgs.callPackage ./examples/moving_poiseulle_particles_2d_hlbm_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; inherit kel; @@ -202,8 +207,6 @@ in rec { examples.poiseulle_particles_2d_ibm_gpu ]; }; - - }; default = release.examples; diff --git a/examples/moving_poiseulle_particles_2d_hlbm_gpu/.nix/derivation.nix b/examples/moving_poiseulle_particles_2d_hlbm_gpu/.nix/derivation.nix new file mode 100644 index 0000000..efc1125 --- /dev/null +++ b/examples/moving_poiseulle_particles_2d_hlbm_gpu/.nix/derivation.nix @@ -0,0 +1,41 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, python3 +, pname +, version +, adaptive-cpp +, kel +}: + +stdenv.mkDerivation { + pname = pname + "-examples-" + "moving_poiseulle_particles_2d_hlbm_gpu"; + inherit version; + src = ./..; + + nativeBuildInputs = [ + scons + clang-tools + python3 + ]; + + buildInputs = [ + forstio.core + forstio.async + forstio.codec + forstio.codec-unit + forstio.io + forstio.remote + forstio.remote-filesystem + forstio.codec-json + adaptive-cpp + kel.lbm.core + kel.lbm.sycl + ]; + + preferLocalBuild = true; + + outputs = [ "out" "dev" ]; +} diff --git a/examples/moving_poiseulle_particles_2d_hlbm_gpu/SConscript b/examples/moving_poiseulle_particles_2d_hlbm_gpu/SConscript new file mode 100644 index 0000000..0953a04 --- /dev/null +++ b/examples/moving_poiseulle_particles_2d_hlbm_gpu/SConscript @@ -0,0 +1,34 @@ +#!/bin/false + +import os +import os.path +import glob + + +Import('env') + +dir_path = Dir('.').abspath + +# Environment for base library +examples_env = env.Clone(); +examples_env['CXX'] = 'syclcc-clang'; +examples_env['CXXFLAGS'] += ['-O3']; + +examples_env.sources = sorted(glob.glob(dir_path + "/*.cpp")) +examples_env.headers = sorted(glob.glob(dir_path + "/*.hpp")) + +env.sources += examples_env.sources; +env.headers += examples_env.headers; + +# Cavity2D +examples_objects = []; +examples_env.add_source_files(examples_objects, ['sim.cpp'], shared=False); +examples_env.poiseulle_2d_gpu = examples_env.Program('#bin/moving_poiseulle_particles_2d_hlbm_gpu', [examples_objects]); + +# Set Alias +env.examples = [ + examples_env.poiseulle_2d_gpu +]; +env.Alias('examples', env.examples); +env.targets += ['examples']; +env.Install('$prefix/bin/', env.examples); diff --git a/examples/moving_poiseulle_particles_2d_hlbm_gpu/SConstruct b/examples/moving_poiseulle_particles_2d_hlbm_gpu/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/moving_poiseulle_particles_2d_hlbm_gpu/SConstruct @@ -0,0 +1,81 @@ +#!/usr/bin/env python3 + +import sys +import os +import os.path +import glob +import re + + +if sys.version_info < (3,): + def isbasestring(s): + return isinstance(s,basestring) +else: + def isbasestring(s): + return isinstance(s, (str,bytes)) + +def add_kel_source_files(self, sources, filetype, lib_env=None, shared=False, target_post=""): + + if isbasestring(filetype): + dir_path = self.Dir('.').abspath + filetype = sorted(glob.glob(dir_path+"/"+filetype)) + + for path in filetype: + target_name = re.sub( r'(.*?)(\.cpp|\.c\+\+)', r'\1' + target_post, path ) + if shared: + target_name+='.os' + sources.append( self.SharedObject( target=target_name, source=path ) ) + else: + target_name+='.o' + sources.append( self.StaticObject( target=target_name, source=path ) ) + pass + +def isAbsolutePath(key, dirname, env): + assert os.path.isabs(dirname), "%r must have absolute path syntax" % (key,) + +env_vars = Variables( + args=ARGUMENTS +) + +env_vars.Add('prefix', + help='Installation target location of build results and headers', + default='/usr/local/', + validator=isAbsolutePath +) + +env_vars.Add('build_examples', + help='If examples should be built', + default="true" +) + +env=Environment(ENV=os.environ, variables=env_vars, CPPPATH=[], + CPPDEFINES=['SAW_UNIX'], + CXXFLAGS=[ + '-std=c++20', + '-g', + '-Wall', + '-Wextra' + ], + LIBS=[ + 'forstio-core', + 'forstio-async', + 'forstio-io' + ] +); +env.__class__.add_source_files = add_kel_source_files +env.Tool('compilation_db'); +env.cdb = env.CompilationDatabase('compile_commands.json'); + +env.objects = []; +env.sources = []; +env.headers = []; +env.targets = []; + +Export('env') +SConscript('SConscript') + +env.Alias('cdb', env.cdb); +env.Alias('all', [env.targets]); +env.Default('all'); + +env.Alias('install', '$prefix') diff --git a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp new file mode 100644 index 0000000..1b7ec2f --- /dev/null +++ b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp @@ -0,0 +1,428 @@ +#include +#include +#include + +#include +#include +#include +#include + +namespace kel { +namespace lbm { + +constexpr uint64_t dim_y = 256ul; +constexpr uint64_t dim_x = dim_y * 20ul; + +// constexpr uint64_t particle_amount = 1ul; + +namespace sch { +using namespace saw::schema; + +using InfoChunk = Chunk; + +template +using DfChunk = Chunk, 1u, dim_x, dim_y>; + +template +using ScalarChunk = Chunk, 0u, dim_x, dim_y>; + +template +using VectorChunk = Chunk, 0u, dim_x, dim_y>; + +template +using ChunkStruct = Struct< + Member, + Member, "dfs">, + Member, "dfs_old">, + Member, "particle_N">, + Member, "particle_D"> +>; + +template +using VelChunk = Chunk, 0u, dim_x, dim_y>; + +template +using RhoChunk = Chunk, 0u, dim_x, dim_y>; + +template +using MacroStruct = Struct< + Member, "velocity">, + Member, "density">, + Member, "porosity"> +>; + +//template +//using ParticleArray = Array< +// Particle +//>; +} + +template +saw::error_or setup_initial_conditions( + 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){ + info_f.at(index).set(1u); + }, + {}, + info_f.get_dims(), + {} + ); + // Fluid + iterator::apply( + [&](auto& index){ + info_f.at(index).set(2u); + }, + {}, + info_f.get_dims(), + {{1u,1u}} + ); + + // Inflow + iterator::apply( + [&](auto& index){ + info_f.at(index).set(3u); + }, + {{0u,0u}}, + {{1u,dim_y}}, + {{0u,1u}} + ); + + // Outflow + iterator::apply( + [&](auto& index){ + info_f.at(index).set(4u); + }, + {{dim_x-1u,0u}}, + {{dim_x, dim_y}}, + {{0u,1u}} + ); + // + auto& df_f = fields.template get<"dfs_old">(); + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"velocity">(); + auto& por_f = macros.template get<"porosity">(); + + iterator::apply( + [&](auto& index){ + auto& df = df_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; + } + auto eq = equilibrium(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims(), + {{1u,1u}} + ); + /* Technically the particle group should handle this case + iterator::apply( + [&](auto& index){ + saw::data> middle, ind_vec; + middle.at({{0u}}) = dim_x * 0.25; + middle.at({{1u}}) = dim_y * 0.5; + + ind_vec.at({{0u}}) = index.at({{0u}}).template cast_to(); + ind_vec.at({{1u}}) = index.at({{1u}}).template cast_to(); + + auto dist = middle - ind_vec; + auto dist_2 = saw::math::dot(dist,dist); + if(dist_2.at({}).get() < dim_y*dim_y*0.01){ + porous_f.at(index).at({}) = 0.0; + } + }, + {},// 0-index + df_f.get_dims() + ); + */ + + + + return saw::make_void(); +} + +template +saw::error_or step( + saw::data>,encode::Sycl>& fields, + saw::data>,encode::Sycl>& macros, + saw::data>& particles + saw::data t_i, + device& dev +){ + auto& q = dev.get_handle(); + auto& info_f = fields.template get<"info">(); + auto& porous_f = macros.template get<"porosity">(); + + // auto coll_ev = + q.submit([&](acpp::sycl::handler& h){ + component> collision{0.65}; + component> bb; + component,encode::Sycl> abb; + + saw::data> rho_b; + rho_b.at({}) = 1.0; + saw::data> vel_b; + vel_b.at({{0u}}) = 0.015; + + component> equi{rho_b,vel_b}; + + component,encode::Sycl> flow_in{ + [&](){ + uint64_t target_t_i = 64u; + if(t_i.get() < target_t_i){ + return 1.0 + (0.0002 / target_t_i) * t_i.get(); + } + return 1.0002; + }() + }; + 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){ + index.at({{i}}).set(idx[i]); + } + + auto info = info_f.at(index); + + switch(info.get()){ + case 0u: + break; + case 1u: + bb.apply(fields,index,t_i); + break; + case 2u: + collision.apply(fields,macros,index,t_i); + break; + case 3u: + flow_in.apply(fields,index,t_i); + // equi.apply(fields,index,t_i); + collision.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + // equi.apply(fields,index,t_i); + collision.apply(fields,macros,index,t_i); + break; + default: + break; + } + }); + }).wait(); + + + // Step + /* + q.submit([&](acpp::sycl::handler& h){ + // h.depends_on(collision_ev); + }).wait(); + */ + + return saw::make_void(); +} +} +} + +template +saw::error_or lbm_main(int argc, char** argv){ + using namespace kel::lbm; + + using dfi = df_info; + + auto eo_lbm_dir = output_directory(); + if(eo_lbm_dir.is_error()){ + return std::move(eo_lbm_dir.get_error()); + } + auto& lbm_dir = eo_lbm_dir.get_value(); + + auto out_dir = lbm_dir / "moving_poiseulle_particles_2d_hlbm_gpu"; + + { + std::error_code ec; + std::filesystem::create_directories(out_dir,ec); + if(ec != std::errc{}){ + return saw::make_error("Could not create output directory"); + } + } + + converter conv { + // delta_x + {{1.0}}, + // delta_t + {{1.0}} + }; + + print_lbm_meta(conv,{0.1},{1e-4},{0.4 * dim_y}); + + // saw::data> meta{{dim_x,dim_y}}; + auto lbm_data_ptr = saw::heap>>(); + auto lbm_macro_data_ptr = saw::heap>>(); + + std::cout<<"Estimated Bytes: "<,sch::MacroStruct>().get()<(*lbm_data_ptr,*lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_vtk_file(out_dir,"initial_state",0u,*lbm_data_ptr); + if(eov.is_error()){ + return eov; + } + } + + saw::data, encode::Sycl> lbm_sycl_data{sycl_q}; + saw::data, encode::Sycl> lbm_sycl_macro_data{sycl_q}; + sycl_q.wait(); + + { + auto eov = dev.copy_to_device(*lbm_data_ptr,lbm_sycl_data); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = dev.copy_to_device(*lbm_macro_data_ptr,lbm_sycl_macro_data); + if(eov.is_error()){ + return eov; + } + } + sycl_q.wait(); + auto lsd_view = make_view(lbm_sycl_data); + auto lsdm_view = make_view(lbm_sycl_macro_data); + saw::data time_steps{16u*4096ul}; + auto& info_f = lsd_view.template get<"info">(); + + for(saw::data i{0u}; i < time_steps and krun; ++i){ + // BC + Collision + { + auto eov = step(lsd_view,lsdm_view,i,dev); + if(eov.is_error()){ + return eov; + } + } + sycl_q.wait(); + if(i.get() % 32u == 0u){ + { + auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_csv_file(out_dir,"m",i.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + } + // Stream + sycl_q.submit([&](acpp::sycl::handler& h){ + component> stream; + + 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]); + } + + auto info = info_f.at(index); + + if(info.get() > 0u){ + stream.apply(lsd_view,index,i); + } + }); + }).wait(); + wait.poll(); + if(print_status){ + std::cout<<"Status: "<().get() * 100 / time_steps.get())<<"%"<(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + std::cerr<<"[Error] "< 0u){ + std::cerr<<" - "< saw::error_or setup_initial_conditions( saw::data>& fields, saw::data>& macros, - saw::data>& particles + saw::data>& p_group ){ auto& info_f = fields.template get<"info">(); // Set everything as walls @@ -207,8 +207,6 @@ saw::error_or step( auto& p_rot = pirb.template get<"rotation">(); auto& p_acc = pirb.template get<"acceleration">(); - - saw::data> delta_t; delta_t.at({}).set(1.0f); verlet_step_lambda(p,delta_t); @@ -218,7 +216,7 @@ saw::error_or step( // auto coll_ev = q.submit([&](acpp::sycl::handler& h){ // Need nicer things to handle the flow. I see improvement here - component> collision{0.8}; + component> collision{0.8}; component> bb; h.parallel_for(acpp::sycl::range{dim_x,dim_y}, [=](acpp::sycl::id idx){ @@ -236,7 +234,7 @@ saw::error_or step( bb.apply(fields,index,t_i); break; case 2u: - // collision.apply(fields,macros,index,t_i); + collision.apply(fields,macros,index,t_i); break; default: break; @@ -268,7 +266,7 @@ saw::error_or lbm_main(int argc, char** argv){ } auto& lbm_env = eo_lbm_env.get_value(); - auto out_dir = lbm_env.data_dir / "settling_cubes_2d_ibm_gpu"; + auto out_dir = lbm_env.data_dir / "settling_spheres_2d_hlbm_gpu"; { std::error_code ec; -- cgit v1.2.3 From 3d021acc9f5283c116d0ec0b5a46d2428e99382c Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Thu, 4 Jun 2026 22:21:45 +0200 Subject: Dangling changes --- default.nix | 5 + .../moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp | 5 +- examples/poiseulle_particles_2d_gpu/common.hpp | 62 +++++ examples/poiseulle_particles_2d_gpu/init.hpp | 122 ++++++++++ examples/poiseulle_particles_2d_gpu/sim.cpp | 259 +-------------------- examples/poiseulle_particles_2d_gpu/step.hpp | 90 +++++++ 6 files changed, 293 insertions(+), 250 deletions(-) create mode 100644 examples/poiseulle_particles_2d_gpu/common.hpp create mode 100644 examples/poiseulle_particles_2d_gpu/init.hpp create mode 100644 examples/poiseulle_particles_2d_gpu/step.hpp diff --git a/default.nix b/default.nix index 99e306a..53c11d3 100644 --- a/default.nix +++ b/default.nix @@ -152,6 +152,11 @@ in rec { 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; + }; + poiseulle_3d = pkgs.callPackage ./examples/poiseulle_3d/.nix/derivation.nix { inherit pname version stdenv forstio adaptive-cpp; inherit kel; diff --git a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp index 1b7ec2f..ce93d3e 100644 --- a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp +++ b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp @@ -55,6 +55,9 @@ using MacroStruct = Struct< //using ParticleArray = Array< // Particle //>; + +template +using ParticleSpheroidGroup = ParticleGroup>; } template @@ -171,7 +174,7 @@ template saw::error_or step( saw::data>,encode::Sycl>& fields, saw::data>,encode::Sycl>& macros, - saw::data>& particles + saw::data>,encode::Sycl>& particles saw::data t_i, device& dev ){ diff --git a/examples/poiseulle_particles_2d_gpu/common.hpp b/examples/poiseulle_particles_2d_gpu/common.hpp new file mode 100644 index 0000000..a69a2cf --- /dev/null +++ b/examples/poiseulle_particles_2d_gpu/common.hpp @@ -0,0 +1,62 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace kel { +namespace lbm { + +constexpr uint64_t dim_y = 256ul; +constexpr uint64_t dim_x = dim_y * 20ul; + +constexpr uint64_t particle_amount = 1ul; + +namespace sch { +using namespace saw::schema; + +using InfoChunk = Chunk; + +template +using DfChunk = Chunk, 1u, dim_x, dim_y>; + +template +using ScalarChunk = Chunk, 0u, dim_x, dim_y>; + +template +using VectorChunk = Chunk, 0u, dim_x, dim_y>; + +template +using ChunkStruct = Struct< + Member, + Member, "dfs">, + Member, "dfs_old">, + Member, "particle_N">, + Member, "particle_D"> +>; + +template +using VelChunk = Chunk, 0u, dim_x, dim_y>; + +template +using RhoChunk = Chunk, 0u, dim_x, dim_y>; + +template +using MacroStruct = Struct< + Member, "velocity">, + Member, "density">, + Member, "porosity"> +>; + +template +using ParticleSpheroidGroup = ParticleGroup>; +} + +} +} diff --git a/examples/poiseulle_particles_2d_gpu/init.hpp b/examples/poiseulle_particles_2d_gpu/init.hpp new file mode 100644 index 0000000..70d59fc --- /dev/null +++ b/examples/poiseulle_particles_2d_gpu/init.hpp @@ -0,0 +1,122 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +template +saw::error_or setup_initial_conditions( + 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){ + info_f.at(index).set(1u); + }, + {}, + info_f.get_dims(), + {} + ); + // Fluid + iterator::apply( + [&](auto& index){ + info_f.at(index).set(2u); + }, + {}, + info_f.get_dims(), + {{1u,1u}} + ); + + // Inflow + iterator::apply( + [&](auto& index){ + info_f.at(index).set(3u); + }, + {{0u,0u}}, + {{1u,dim_y}}, + {{0u,1u}} + ); + + // Outflow + iterator::apply( + [&](auto& index){ + info_f.at(index).set(4u); + }, + {{dim_x-1u,0u}}, + {{dim_x, dim_y}}, + {{0u,1u}} + ); + // + auto& df_f = fields.template get<"dfs_old">(); + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"velocity">(); + auto& por_f = macros.template get<"porosity">(); + + iterator::apply( + [&](auto& index){ + auto& df = df_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; + } + auto eq = equilibrium(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims(), + {{1u,1u}} + ); + + iterator::apply( + [&](auto& index){ + saw::data> middle, ind_vec; + middle.at({{0u}}) = dim_x * 0.25; + middle.at({{1u}}) = dim_y * 0.5; + + ind_vec.at({{0u}}) = index.at({{0u}}).template cast_to(); + ind_vec.at({{1u}}) = index.at({{1u}}).template cast_to(); + + auto dist = middle - ind_vec; + auto dist_2 = saw::math::dot(dist,dist); + if(dist_2.at({}).get() < dim_y*dim_y*0.01){ + porous_f.at(index).at({}) = 0.0; + } + }, + {},// 0-index + df_f.get_dims() + ); + + { + saw::data> dense_p; + dense_p.at({}).set(1); + particles = create_spheroid_particle_group(dense_p, {{16u}}); + } + + return saw::make_void(); +} + +} +} diff --git a/examples/poiseulle_particles_2d_gpu/sim.cpp b/examples/poiseulle_particles_2d_gpu/sim.cpp index 9375078..756c5b4 100644 --- a/examples/poiseulle_particles_2d_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_gpu/sim.cpp @@ -1,249 +1,6 @@ -#include -#include -#include - -#include -#include -#include -#include - -namespace kel { -namespace lbm { - -constexpr uint64_t dim_y = 256ul; -constexpr uint64_t dim_x = dim_y * 20ul; - -constexpr uint64_t particle_amount = 1ul; - -namespace sch { -using namespace saw::schema; - -using InfoChunk = Chunk; - -template -using DfChunk = Chunk, 1u, dim_x, dim_y>; - -template -using ScalarChunk = Chunk, 0u, dim_x, dim_y>; - -template -using VectorChunk = Chunk, 0u, dim_x, dim_y>; - -template -using ChunkStruct = Struct< - Member, - Member, "dfs">, - Member, "dfs_old">, - Member, "particle_N">, - Member, "particle_D"> ->; - -template -using VelChunk = Chunk, 0u, dim_x, dim_y>; - -template -using RhoChunk = Chunk, 0u, dim_x, dim_y>; - -template -using MacroStruct = Struct< - Member, "velocity">, - Member, "density">, - Member, "porosity"> ->; - -//template -//using ParticleArray = Array< -// Particle -//>; -} - -template -saw::error_or setup_initial_conditions( - saw::data>& fields, - saw::data>& macros -){ - auto& info_f = fields.template get<"info">(); - auto& porous_f = macros.template get<"porosity">(); - // Set everything as walls - iterator::apply( - [&](auto& index){ - info_f.at(index).set(1u); - }, - {}, - info_f.get_dims(), - {} - ); - // Fluid - iterator::apply( - [&](auto& index){ - info_f.at(index).set(2u); - }, - {}, - info_f.get_dims(), - {{1u,1u}} - ); - - // Inflow - iterator::apply( - [&](auto& index){ - info_f.at(index).set(3u); - }, - {{0u,0u}}, - {{1u,dim_y}}, - {{0u,1u}} - ); - - // Outflow - iterator::apply( - [&](auto& index){ - info_f.at(index).set(4u); - }, - {{dim_x-1u,0u}}, - {{dim_x, dim_y}}, - {{0u,1u}} - ); - // - auto& df_f = fields.template get<"dfs_old">(); - auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); - auto& por_f = macros.template get<"porosity">(); - - iterator::apply( - [&](auto& index){ - auto& df = df_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; - } - auto eq = equilibrium(rho,vel); - - df = eq; - }, - {},// 0-index - df_f.get_dims(), - {{1u,1u}} - ); - - iterator::apply( - [&](auto& index){ - saw::data> middle, ind_vec; - middle.at({{0u}}) = dim_x * 0.25; - middle.at({{1u}}) = dim_y * 0.5; - - ind_vec.at({{0u}}) = index.at({{0u}}).template cast_to(); - ind_vec.at({{1u}}) = index.at({{1u}}).template cast_to(); - - auto dist = middle - ind_vec; - auto dist_2 = saw::math::dot(dist,dist); - if(dist_2.at({}).get() < dim_y*dim_y*0.01){ - porous_f.at(index).at({}) = 0.0; - } - }, - {},// 0-index - df_f.get_dims() - ); - - return saw::make_void(); -} - -template -saw::error_or step( - saw::data>,encode::Sycl>& fields, - saw::data>,encode::Sycl>& macros, - saw::data t_i, - device& dev -){ - auto& q = dev.get_handle(); - auto& info_f = fields.template get<"info">(); - auto& porous_f = macros.template get<"porosity">(); - - // auto coll_ev = - q.submit([&](acpp::sycl::handler& h){ - component> collision{0.65}; - component> bb; - component,encode::Sycl> abb; - - saw::data> rho_b; - rho_b.at({}) = 1.0; - saw::data> vel_b; - vel_b.at({{0u}}) = 0.015; - - component> equi{rho_b,vel_b}; - - component,encode::Sycl> flow_in{ - [&](){ - uint64_t target_t_i = 64u; - if(t_i.get() < target_t_i){ - return 1.0 + (0.0002 / target_t_i) * t_i.get(); - } - return 1.0002; - }() - }; - 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){ - index.at({{i}}).set(idx[i]); - } - - auto info = info_f.at(index); - - switch(info.get()){ - case 0u: - break; - case 1u: - bb.apply(fields,index,t_i); - break; - case 2u: - collision.apply(fields,macros,index,t_i); - break; - case 3u: - flow_in.apply(fields,index,t_i); - // equi.apply(fields,index,t_i); - collision.apply(fields,macros,index,t_i); - break; - case 4u: - flow_out.apply(fields,index,t_i); - // equi.apply(fields,index,t_i); - collision.apply(fields,macros,index,t_i); - break; - default: - break; - } - }); - }).wait(); - - - // Step - /* - q.submit([&](acpp::sycl::handler& h){ - // h.depends_on(collision_ev); - }).wait(); - */ - - return saw::make_void(); -} -} -} +#include "common.hpp" +#include "init.hpp" +#include "step.hpp" template saw::error_or lbm_main(int argc, char** argv){ @@ -257,7 +14,7 @@ saw::error_or lbm_main(int argc, char** argv){ } auto& lbm_dir = eo_lbm_dir.get_value(); - auto out_dir = lbm_dir / "poiseulle_particles_2d_hlbm_gpu"; + auto out_dir = lbm_dir / "poiseulle_particles_2d_gpu"; { std::error_code ec; @@ -279,6 +36,7 @@ saw::error_or lbm_main(int argc, char** argv){ // saw::data> meta{{dim_x,dim_y}}; auto lbm_data_ptr = saw::heap>>(); auto lbm_macro_data_ptr = saw::heap>>(); + auto lbm_particle_data_ptr = saw::heap>>(); std::cout<<"Estimated Bytes: "<,sch::MacroStruct>().get()< lbm_main(int argc, char** argv){ sycl_q.wait(); { - auto eov = setup_initial_conditions(*lbm_data_ptr,*lbm_macro_data_ptr); + auto eov = setup_initial_conditions(*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_data_ptr); if(eov.is_error()){ return eov; } @@ -318,6 +76,7 @@ saw::error_or lbm_main(int argc, char** argv){ saw::data, encode::Sycl> lbm_sycl_data{sycl_q}; saw::data, encode::Sycl> lbm_sycl_macro_data{sycl_q}; + saw::data, encode::Sycl> lbm_sycl_particle_data{sycl_q}; sycl_q.wait(); { @@ -335,13 +94,15 @@ saw::error_or lbm_main(int argc, char** argv){ sycl_q.wait(); auto lsd_view = make_view(lbm_sycl_data); auto lsdm_view = make_view(lbm_sycl_macro_data); + auto lsdp_view = make_view(lbm_sycl_particle_data); + saw::data time_steps{16u*4096ul}; auto& info_f = lsd_view.template get<"info">(); for(saw::data i{0u}; i < time_steps and krun; ++i){ // BC + Collision { - auto eov = step(lsd_view,lsdm_view,i,dev); + auto eov = step(lsd_view,lsdm_view,lsdp_view,i,dev); if(eov.is_error()){ return eov; } diff --git a/examples/poiseulle_particles_2d_gpu/step.hpp b/examples/poiseulle_particles_2d_gpu/step.hpp new file mode 100644 index 0000000..a4e44b4 --- /dev/null +++ b/examples/poiseulle_particles_2d_gpu/step.hpp @@ -0,0 +1,90 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +template +saw::error_or step( + saw::data>,encode::Sycl>& fields, + saw::data>,encode::Sycl>& macros, + saw::data>,encode::Sycl>& particles, + saw::data t_i, + device& dev +){ + auto& q = dev.get_handle(); + auto& info_f = fields.template get<"info">(); + auto& porous_f = macros.template get<"porosity">(); + + // auto coll_ev = + q.submit([&](acpp::sycl::handler& h){ + component> collision{0.65}; + component> bb; + component,encode::Sycl> abb; + + saw::data> rho_b; + rho_b.at({}) = 1.0; + saw::data> vel_b; + vel_b.at({{0u}}) = 0.015; + + component> equi{rho_b,vel_b}; + + component,encode::Sycl> flow_in{ + [&](){ + uint64_t target_t_i = 64u; + if(t_i.get() < target_t_i){ + return 1.0 + (0.0002 / target_t_i) * t_i.get(); + } + return 1.0002; + }() + }; + 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){ + index.at({{i}}).set(idx[i]); + } + + auto info = info_f.at(index); + + switch(info.get()){ + case 0u: + break; + case 1u: + bb.apply(fields,index,t_i); + break; + case 2u: + collision.apply(fields,macros,index,t_i); + break; + case 3u: + flow_in.apply(fields,index,t_i); + // equi.apply(fields,index,t_i); + collision.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + // equi.apply(fields,index,t_i); + collision.apply(fields,macros,index,t_i); + break; + default: + break; + } + }); + }).wait(); + + + // Step + /* + q.submit([&](acpp::sycl::handler& h){ + // h.depends_on(collision_ev); + }).wait(); + */ + + return saw::make_void(); +} + +} +} -- cgit v1.2.3 From 932fbf86d60df48623ad5fbc9d60e572bb68ef12 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Mon, 8 Jun 2026 20:15:20 +0200 Subject: Dangling changes --- examples/poiseulle_particles_2d_gpu/sim.cpp | 2 +- lib/core/c++/particle/particle.hpp | 3 --- lib/core/c++/simulation/common.hpp | 7 +++++++ lib/core/c++/simulation/hlbm.hpp | 12 ++++++++++++ 4 files changed, 20 insertions(+), 4 deletions(-) create mode 100644 lib/core/c++/simulation/common.hpp create mode 100644 lib/core/c++/simulation/hlbm.hpp diff --git a/examples/poiseulle_particles_2d_gpu/sim.cpp b/examples/poiseulle_particles_2d_gpu/sim.cpp index 756c5b4..3de3cfb 100644 --- a/examples/poiseulle_particles_2d_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_gpu/sim.cpp @@ -14,7 +14,7 @@ saw::error_or lbm_main(int argc, char** argv){ } auto& lbm_dir = eo_lbm_dir.get_value(); - auto out_dir = lbm_dir / "poiseulle_particles_2d_gpu"; + auto out_dir = lbm_dir / "poiseulle_particles_2d_gpu" / "hlbm"; { std::error_code ec; diff --git a/lib/core/c++/particle/particle.hpp b/lib/core/c++/particle/particle.hpp index fec2eca..1a99dcd 100644 --- a/lib/core/c++/particle/particle.hpp +++ b/lib/core/c++/particle/particle.hpp @@ -62,8 +62,6 @@ using ParticleGroup = Struct< >; } - - template::type radius> saw::data>> create_spheroid_particle_group( saw::data> density_p, @@ -287,7 +285,6 @@ constexpr auto broadphase_collision_check = []( namespace impl { - } } diff --git a/lib/core/c++/simulation/common.hpp b/lib/core/c++/simulation/common.hpp new file mode 100644 index 0000000..f675a99 --- /dev/null +++ b/lib/core/c++/simulation/common.hpp @@ -0,0 +1,7 @@ +#pragma once + +namespace kel { +namespace lbm { + +} +} diff --git a/lib/core/c++/simulation/hlbm.hpp b/lib/core/c++/simulation/hlbm.hpp new file mode 100644 index 0000000..93df9a7 --- /dev/null +++ b/lib/core/c++/simulation/hlbm.hpp @@ -0,0 +1,12 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +class simulation {}; + + +} +} -- cgit v1.2.3