From e85b713cf31697a3309e12f30ba5759fee1cd3cc Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Fri, 20 Mar 2026 16:08:25 +0100 Subject: Trying out n linear interpolation --- lib/core/c++/math/n_linear.hpp | 44 +++++++++++++++++++++++++++++------------- lib/core/tests/math.cpp | 20 +++++++++++++++++++ 2 files changed, 51 insertions(+), 13 deletions(-) create mode 100644 lib/core/tests/math.cpp diff --git a/lib/core/c++/math/n_linear.hpp b/lib/core/c++/math/n_linear.hpp index 62906cd..4f47965 100644 --- a/lib/core/c++/math/n_linear.hpp +++ b/lib/core/c++/math/n_linear.hpp @@ -6,24 +6,42 @@ namespace kel { namespace lbm { template -void n_linear_interpolate(const saw::data& field, const saw::data>& pos){ +auto n_linear_interpolate(const saw::data& field, const saw::data>& pos){ + // Pos auto pos_bound = pos; + + // Dimensions auto meta = field.dims(); - saw::data ind; + + // Lower Index + saw::data> ind; + for(saw::data i{0u}; i < saw::data{D}; ++i){ - pos_bound.at(i).set( - std::max( - std::min( - meta.at(i).template cast_to().get(), - pos.at(i).get() + 1.5 - ), - 1, - ) - -1 - ); - ind.at(i) = pos_bound.at(i).template cast_to(); + // Native Positive i + auto npos_i = pos.at({i}).get(); + + { + // Ok I want to stay in bounds + npos_i = std::min(npos_i,meta.at(i).get()-1.0); + npos_i = std::max(npos_i,1.0); + } + + // Native Index i + auto nind_i = static_cast(std::floor(npos_i))-1ul; + + // Set index to i + ind.at(i).set(nind_i); } + saw::data> pos_frac; + for(saw::data i{0u}; i < saw::data{D}; ++i){ + pos_frac.at({i}) = pos_bound.at({i}) - ind.at(i).template cast_to(); + } + + // Base value + auto res = field.at(ind); + + constexpr uint64_t d_corners = 1ull << D; } } } diff --git a/lib/core/tests/math.cpp b/lib/core/tests/math.cpp new file mode 100644 index 0000000..738d637 --- /dev/null +++ b/lib/core/tests/math.cpp @@ -0,0 +1,20 @@ +#include + +#include "../c++/math/n_linear.hpp" + +namespace { +namespace sch { +using namespace saw::schema; +} + +SAW_TEST("Math 2-Linear"){ + using namespace kel; + + saw::data,32u,32u>> field; + + saw::data> pos; + + SAW_EXPECT(true, "Default true check"); +} + +} -- cgit v1.2.3 From 16c1198f8fb401c3a98d927053fb2d29c2ce5f91 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Fri, 20 Mar 2026 18:17:55 +0100 Subject: Too much time spent on downcasts, but well. Here I am. --- lib/core/c++/math/n_linear.hpp | 76 ++++++++++++++++++++++++++++++++++++++++-- lib/core/tests/math.cpp | 26 +++++++++++++-- 2 files changed, 97 insertions(+), 5 deletions(-) diff --git a/lib/core/c++/math/n_linear.hpp b/lib/core/c++/math/n_linear.hpp index 4f47965..bb7653c 100644 --- a/lib/core/c++/math/n_linear.hpp +++ b/lib/core/c++/math/n_linear.hpp @@ -1,12 +1,51 @@ #pragma once #include "../common.hpp" +#include "../iterator.hpp" namespace kel { namespace lbm { +namespace impl { +template +struct n_linear_interpolate_helper final { + template + auto apply(const saw::data& field, const saw::data>& pos){ + return pos; + } +}; +} + +template +saw::data,sch::Vector>> position_to_index_and_fraction(const saw::data>& pos){ + saw::data,sch::Vector>> sep; + + auto& ind = sep.template get<0u>(); + auto& frac = sep.template get<1u>(); + + auto pos_cpy = pos; + // Guarantee that the pos is at least 0 + for(uint64_t i = 0u; i < D; ++i){ + pos_cpy.at({{i}}).set(std::max(pos.at({{i}}).get(), static_cast::type>(0))); + } + + // Now we can cast to uint64_t + for(uint64_t i = 0u; i < D; ++i){ + ind.at({{i}}) = pos_cpy.at({{i}}).template cast_to(); + } + + frac = pos - ind.template cast_to(); + + return sep; +} + +template +auto floor_index_from_position(const saw::data>& pos){ + return position_to_index_and_fraction(pos).template get<0u>(); +} template -auto n_linear_interpolate(const saw::data& field, const saw::data>& pos){ +auto n_linear_interpolate( + const saw::data& field, const saw::data>& pos){ // Pos auto pos_bound = pos; @@ -39,9 +78,40 @@ auto n_linear_interpolate(const saw::data& field, const saw::data res; - constexpr uint64_t d_corners = 1ull << D; + // constexpr uint64_t d_corners = 1ul << D; + + saw::data> ones_ind; + for(saw::data i{0u}; i < saw::data{D}; ++i){ + ones_ind.at({i}).set(1u); + } + + iterator::apply([&](auto ind){ + // Iterates over (0,0,0) to (1,1,1) + saw::data weight{1.0}; + + for(saw::data d{0u}; d < saw::data{D}; ++d){ + + saw::data t = pos_frac.at({d}); + + if(ind.at(d).get() == 0u){ + weight = weight * (saw::data{1} - t); + }else{ + weight = weight * t; + } + } + + }, {}, ones_ind); +} + +template +saw::data> bilinear_interpolate(const saw::data& field, const saw::data>& pos){ + saw::data> res; + + { + } + return {}; } } } diff --git a/lib/core/tests/math.cpp b/lib/core/tests/math.cpp index 738d637..7c9dc14 100644 --- a/lib/core/tests/math.cpp +++ b/lib/core/tests/math.cpp @@ -1,5 +1,6 @@ #include +#include #include "../c++/math/n_linear.hpp" namespace { @@ -7,13 +8,34 @@ namespace sch { using namespace saw::schema; } -SAW_TEST("Math 2-Linear"){ +SAW_TEST("Math 1-Linear"){ using namespace kel; + + saw::data,2u>> field; + { + field.at({{0u}}).at({{0u}}).set(-1.0); + field.at({{1u}}).at({{0u}}).set(2.0); + } + saw::data> pos; + pos.at({{0u}}).set(0.3); +} - saw::data,32u,32u>> field; +SAW_TEST("Math/Floor Index from Pos"){ + using namespace kel; saw::data> pos; + { + pos.at({{0u}}) = 43.999; + pos.at({{1u}}) = -50.0; + } + + auto ind = lbm::floor_index_from_position(pos); + for(uint64_t i = 0u; i < 2u; ++i){ + std::cout< Date: Fri, 20 Mar 2026 18:23:39 +0100 Subject: Fraction and downcast looks ok. Need upper bound now --- lib/core/c++/math/n_linear.hpp | 2 +- lib/core/tests/math.cpp | 8 +++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/lib/core/c++/math/n_linear.hpp b/lib/core/c++/math/n_linear.hpp index bb7653c..ac02da8 100644 --- a/lib/core/c++/math/n_linear.hpp +++ b/lib/core/c++/math/n_linear.hpp @@ -33,7 +33,7 @@ saw::data,sch::Vector>> position_to_i ind.at({{i}}) = pos_cpy.at({{i}}).template cast_to(); } - frac = pos - ind.template cast_to(); + frac = pos_cpy - ind.template cast_to(); return sep; } diff --git a/lib/core/tests/math.cpp b/lib/core/tests/math.cpp index 7c9dc14..970004b 100644 --- a/lib/core/tests/math.cpp +++ b/lib/core/tests/math.cpp @@ -30,11 +30,17 @@ SAW_TEST("Math/Floor Index from Pos"){ pos.at({{1u}}) = -50.0; } - auto ind = lbm::floor_index_from_position(pos); + auto ind_frac = lbm::position_to_index_and_fraction(pos); + auto& ind = ind_frac.template get<0u>(); for(uint64_t i = 0u; i < 2u; ++i){ std::cout<(); + for(uint64_t i = 0u; i < 2u; ++i){ + std::cout< Date: Fri, 20 Mar 2026 18:33:31 +0100 Subject: Wanted to finish n linearity today, but I guess I should be happy with this --- lib/core/c++/math/n_linear.hpp | 19 +++++++++++++++++++ lib/core/tests/math.cpp | 33 ++++++++++++++++++++++++++++++++- 2 files changed, 51 insertions(+), 1 deletion(-) diff --git a/lib/core/c++/math/n_linear.hpp b/lib/core/c++/math/n_linear.hpp index ac02da8..8fb0600 100644 --- a/lib/core/c++/math/n_linear.hpp +++ b/lib/core/c++/math/n_linear.hpp @@ -43,6 +43,25 @@ auto floor_index_from_position(const saw::data>& pos){ return position_to_index_and_fraction(pos).template get<0u>(); } +template +saw::data,sch::Vector>> position_to_index_and_fraction_bounded( + const saw::data>& pos, + const saw::data>& bound) +{ + auto infr = position_to_index_and_fraction(pos); + auto& ind = infr.template get<0u>(); + auto& fra = infr.template get<1u>(); + for(uint64_t i = 0u; i < D; ++i){ + // If index is higher than bound. Set to bound and reset fraction + if((ind.at({{i}}).get()+1u) >= bound.at({{i}}).get()){ + ind.at({{i}}).set(bound.at({{i}}).get()-1u); + fra.at({{i}}) = {}; + } + } + return infr; +} + + template auto n_linear_interpolate( const saw::data& field, const saw::data>& pos){ diff --git a/lib/core/tests/math.cpp b/lib/core/tests/math.cpp index 970004b..d456ce8 100644 --- a/lib/core/tests/math.cpp +++ b/lib/core/tests/math.cpp @@ -20,7 +20,7 @@ SAW_TEST("Math 1-Linear"){ pos.at({{0u}}).set(0.3); } -SAW_TEST("Math/Floor Index from Pos"){ +SAW_TEST("Math/Floor Index and Fraction from Position"){ using namespace kel; saw::data> pos; @@ -45,4 +45,35 @@ SAW_TEST("Math/Floor Index from Pos"){ SAW_EXPECT(true, "Default true check"); } +SAW_TEST("Math/Floor Index and Fraction from Position and Bounded"){ + using namespace kel; + + saw::data> pos; + + { + pos.at({{0u}}) = 43.999; + pos.at({{1u}}) = -50.0; + } + + saw::data> bound; + { + bound.at({{0u}}) = 32u; + bound.at({{1u}}) = 16u; + } + + auto ind_frac = lbm::position_to_index_and_fraction_bounded(pos,bound); + auto& ind = ind_frac.template get<0u>(); + for(uint64_t i = 0u; i < 2u; ++i){ + std::cout<(); + for(uint64_t i = 0u; i < 2u; ++i){ + std::cout< Date: Sun, 22 Mar 2026 12:55:23 +0100 Subject: Starting settling cube example --- examples/settling_cubes_2d_ibm/.nix/derivation.nix | 41 ++ examples/settling_cubes_2d_ibm/SConscript | 34 ++ examples/settling_cubes_2d_ibm/SConstruct | 81 ++++ examples/settling_cubes_2d_ibm/sim.cpp | 454 +++++++++++++++++++++ 4 files changed, 610 insertions(+) create mode 100644 examples/settling_cubes_2d_ibm/.nix/derivation.nix create mode 100644 examples/settling_cubes_2d_ibm/SConscript create mode 100644 examples/settling_cubes_2d_ibm/SConstruct create mode 100644 examples/settling_cubes_2d_ibm/sim.cpp diff --git a/examples/settling_cubes_2d_ibm/.nix/derivation.nix b/examples/settling_cubes_2d_ibm/.nix/derivation.nix new file mode 100644 index 0000000..d7f138b --- /dev/null +++ b/examples/settling_cubes_2d_ibm/.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_cubes_2d_ibm_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_cubes_2d_ibm/SConscript b/examples/settling_cubes_2d_ibm/SConscript new file mode 100644 index 0000000..ae08056 --- /dev/null +++ b/examples/settling_cubes_2d_ibm/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_cubes_2d_ibm/SConstruct b/examples/settling_cubes_2d_ibm/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/settling_cubes_2d_ibm/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_cubes_2d_ibm/sim.cpp b/examples/settling_cubes_2d_ibm/sim.cpp new file mode 100644 index 0000000..411efa7 --- /dev/null +++ b/examples/settling_cubes_2d_ibm/sim.cpp @@ -0,0 +1,454 @@ +#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"> +>; + +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"> +>; + +//template +//using ParticleArray = Array< +// Particle +//>; +} + +template +saw::error_or setup_initial_conditions( + saw::data>& fields, + saw::data>& macros, + saw::data, particle_amount>>& 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}} + ); + + // 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">(); + + 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() + ); + + 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){ + info_f.at(index) = 1u; + } + }, + {},// 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, particle_amount>>& particles, + saw::data t_i, + device& dev +){ + auto& q = dev.get_handle(); + auto& info_f = fields.template get<"info">(); + + { + auto& p = particles.at({{0u}}); + + auto& p_coll = p.template get<"collision">(); + auto& p_rad = p_coll.template get<"radius">(); + } + + + // auto coll_ev = + q.submit([&](acpp::sycl::handler& h){ + // Need nicer things to handle the flow. I see improvement here + 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.0015 / target_t_i) * t_i.get(); + } + return 1.0015; + }() + }; + 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: + abb.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 / "poiseulle_particles_2d_bgk_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, particle_amount>>>(); + + std::cout<<"Estimated Bytes: "<,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, particle_amount>, encode::Sycl> lbm_sycl_particle_data{sycl_q}; + sycl_q.wait(); + + auto lsd_view = make_chunk_struct_view(lbm_sycl_data); + auto lsdm_view = make_chunk_struct_view(lbm_sycl_macro_data); + { + 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(); + 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,*lbm_particle_data_ptr,i,dev); + if(eov.is_error()){ + return eov; + } + } + sycl_q.wait(); + /* + { + { + 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<<" - "< Date: Mon, 23 Mar 2026 09:12:42 +0100 Subject: Adding cube generation --- lib/core/c++/particle/geometry/cube.hpp | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) create mode 100644 lib/core/c++/particle/geometry/cube.hpp diff --git a/lib/core/c++/particle/geometry/cube.hpp b/lib/core/c++/particle/geometry/cube.hpp new file mode 100644 index 0000000..6392de8 --- /dev/null +++ b/lib/core/c++/particle/geometry/cube.hpp @@ -0,0 +1,22 @@ +#pragma once + +#include "../particle.hpp" + +namespace kel { +namespace lbm { +template +class particle_cubic_geometry { +private: +public: + template + saw::data> generate_mask(uint64_t resolution, uint64_t bound_nodes = 0u){ + + saw::data> mask; + + auto& grid = mask.template get<"grid">(); + auto& com = mask.template get<"center_of_mass">(); + + + } +} +} -- cgit v1.2.3 From 0d3c1a0f68792a04129dc8da8e9b9b113f76e3ec Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Mon, 23 Mar 2026 16:41:32 +0100 Subject: Found wrong zou he boundary issue. Preparing for fix --- examples/poiseulle_particles_2d_psm_gpu/sim.cpp | 20 +++++++++----------- lib/core/c++/boundary.hpp | 7 ++++--- lib/core/c++/descriptor.hpp | 18 +++++++++--------- 3 files changed, 22 insertions(+), 23 deletions(-) diff --git a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp index a8aeabb..0d356dd 100644 --- a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp @@ -202,7 +202,7 @@ saw::error_or step( q.submit([&](acpp::sycl::handler& h){ - component> collision{0.65}; + component> collision{0.8}; component> bb; component,encode::Sycl> abb; @@ -243,13 +243,13 @@ saw::error_or step( collision.apply(fields,macros,index,t_i); break; case 3u: - //flow_in.apply(fields,index,t_i); - equi.apply(fields,index,t_i); + 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); + flow_out.apply(fields,index,t_i); + // equi.apply(fields,index,t_i); collision.apply(fields,macros,index,t_i); break; default: @@ -295,12 +295,12 @@ saw::error_or lbm_main(int argc, char** argv){ converter conv { // delta_x - {{1.0}}, + {{1e-3}}, // delta_t - {{1.0}} + {{1e-06}} }; - print_lbm_meta(conv,{0.05},{0.01},{0.4 * dim_y}); + print_lbm_meta(conv,{0.1},{1e-2},{0.4 * dim_y *1e-3}); // saw::data> meta{{dim_x,dim_y}}; auto lbm_data_ptr = saw::heap>>(); @@ -382,8 +382,7 @@ saw::error_or lbm_main(int argc, char** argv){ } } sycl_q.wait(); - /* - { + if(i.get() % 1 == 0u){ { auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr); if(eov.is_error()){ @@ -397,7 +396,6 @@ saw::error_or lbm_main(int argc, char** argv){ } } } - */ /* { { diff --git a/lib/core/c++/boundary.hpp b/lib/core/c++/boundary.hpp index adb473d..bb8eed9 100644 --- a/lib/core/c++/boundary.hpp +++ b/lib/core/c++/boundary.hpp @@ -206,16 +206,17 @@ public: /** * Get the sum of the unknown dfs and precalculate the direction */ - auto sum_unknown_dfs = (rho_setting_ - sum_df) * saw::data{known_dir}; + auto sum_unknown_dfs = (rho_setting_ - sum_df); + auto unknown_dfs_dir = sum_unknown_dfs * saw::data{known_dir}; for(saw::data k{0u}; k < saw::data{Descriptor::Q}; ++k){ auto c_k = dfi::directions[k.get()]; if(c_k[0u]*known_dir < 0){ - sum_unknown_dfs += dfs_old.at({k}) * c_k[0u]; + unknown_dfs_dir += dfs_old.at({k}) * c_k[0u]; } } - auto vel_x = sum_unknown_dfs / rho_setting_; + auto vel_x = unknown_dfs_dir / rho_setting_; static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly"); diff --git a/lib/core/c++/descriptor.hpp b/lib/core/c++/descriptor.hpp index 73f0cce..9f7399a 100644 --- a/lib/core/c++/descriptor.hpp +++ b/lib/core/c++/descriptor.hpp @@ -146,15 +146,15 @@ public: static constexpr uint64_t Q = 9u; static constexpr std::array, Q> directions = {{ - { 0, 0}, - {-1, 0}, - { 1, 0}, - { 0,-1}, - { 0, 1}, - {-1,-1}, - { 1, 1}, - {-1, 1}, - { 1,-1} + { 0, 0}, // 0 + {-1, 0}, // 1 + { 1, 0}, // 2 + { 0,-1}, // 3 + { 0, 1}, // 4 + {-1,-1}, // 5 + { 1, 1}, // 6 + {-1, 1}, // 7 + { 1,-1} // 8 }}; static constexpr std::array::type,Q> weights = { -- cgit v1.2.3 From 51fb97c0fc981f5a1e08f83b12f9d0712c69ead3 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Mon, 23 Mar 2026 16:57:33 +0100 Subject: Fixed boundaries for D2Q9 Pressure Zou He --- lib/core/c++/boundary.hpp | 50 +++++++++++++++++++---------------------------- 1 file changed, 20 insertions(+), 30 deletions(-) diff --git a/lib/core/c++/boundary.hpp b/lib/core/c++/boundary.hpp index bb8eed9..ef468be 100644 --- a/lib/core/c++/boundary.hpp +++ b/lib/core/c++/boundary.hpp @@ -191,43 +191,33 @@ public: auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); auto& dfs_old = dfs_old_f.at(index); - /** - * Sum all known DFs - */ - saw::data sum_df{0}; - for(saw::data k{0u}; k < saw::data{Descriptor::Q}; ++k){ - auto c_k = dfi::directions[k.get()]; - - if(c_k[0u]*known_dir <= 0){ - sum_df += dfs_old.at({k}); - } - } - - /** - * Get the sum of the unknown dfs and precalculate the direction - */ - auto sum_unknown_dfs = (rho_setting_ - sum_df); - auto unknown_dfs_dir = sum_unknown_dfs * saw::data{known_dir}; - for(saw::data k{0u}; k < saw::data{Descriptor::Q}; ++k){ - auto c_k = dfi::directions[k.get()]; - if(c_k[0u]*known_dir < 0){ - unknown_dfs_dir += dfs_old.at({k}) * c_k[0u]; + auto rho_vel_x = [&](){ + if constexpr (Dir) -> saw::data { + auto S = dfs_old.at({0u}) + + dfs_old.at({3u}) + dfs_old.at({4u}) + + (dfs_old.at({1u}) + dfs_old.at({5u}) + dfs_old.at({7u})) * 2u; + return rho_setting_ - S; } - } - - auto vel_x = unknown_dfs_dir / rho_setting_; + else if constexpr (not Dir) { + auto S = dfs_old.at({0u}) + + dfs_old.at({3u}) + dfs_old.at({4u}) + + (dfs_old.at({2u}) + dfs_old.at({6u}) + dfs_old.at({8u})) * 2u; + return S - rho_setting_; + } + return {}; + }(); static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly"); if constexpr (Dir) { - dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data{2.0 / 3.0} * rho_setting_ * vel_x; - dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); - dfs_old.at({8u}) = dfs_old.at({7u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); + dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data{2.0 / 3.0} * rho_vel_x; + dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data{1.0 / 6.0} * rho_vel_x + saw::data{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); + dfs_old.at({8u}) = dfs_old.at({7u}) + saw::data{1.0 / 6.0} * rho_vel_x + saw::data{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); }else if constexpr (not Dir){ - dfs_old.at({1u}) = dfs_old.at({2u}) - saw::data{2.0 / 3.0} * rho_setting_ * vel_x; - dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); - dfs_old.at({7u}) = dfs_old.at({8u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); + dfs_old.at({1u}) = dfs_old.at({2u}) - saw::data{2.0 / 3.0} * rho_vel_x; + dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data{1.0 / 6.0} * rho_vel_x + saw::data{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); + dfs_old.at({7u}) = dfs_old.at({8u}) - saw::data{1.0 / 6.0} * rho_vel_x + saw::data{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); } } }; -- cgit v1.2.3 From 46b49e3b4fa283590d120703f80f892ee1f03ffc Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Mon, 23 Mar 2026 17:06:35 +0100 Subject: Fixing wrong lambda return indicator --- lib/core/c++/boundary.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/core/c++/boundary.hpp b/lib/core/c++/boundary.hpp index ef468be..217013d 100644 --- a/lib/core/c++/boundary.hpp +++ b/lib/core/c++/boundary.hpp @@ -192,8 +192,8 @@ public: auto& dfs_old = dfs_old_f.at(index); - auto rho_vel_x = [&](){ - if constexpr (Dir) -> saw::data { + auto rho_vel_x = [&]() -> saw::data { + if constexpr (Dir){ auto S = dfs_old.at({0u}) + dfs_old.at({3u}) + dfs_old.at({4u}) + (dfs_old.at({1u}) + dfs_old.at({5u}) + dfs_old.at({7u})) * 2u; -- cgit v1.2.3 From 571e79c4d0b72202186fd11314cf268723b1844d Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Mon, 23 Mar 2026 17:33:05 +0100 Subject: Fixed boundary setting --- examples/poiseulle_particles_2d_psm_gpu/sim.cpp | 25 +++++++++++++++++++++++-- lib/core/c++/boundary.hpp | 5 ++--- 2 files changed, 25 insertions(+), 5 deletions(-) diff --git a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp index 0d356dd..7702808 100644 --- a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp @@ -83,7 +83,24 @@ saw::error_or setup_initial_conditions( info_f.get_dims(), {{1u,1u}} ); - + // Corners + /// Inflow + iterator::apply( + [&](auto& index){ + info_f.at(index).set(5u); + }, + {{0u,0u}}, + {{1u,dim_y}} + ); + /// Outflow + iterator::apply( + [&](auto& index){ + info_f.at(index).set(5u); + }, + {{dim_x-1u,0u}}, + {{dim_x, dim_y}} + ); + // Overwrite with // Inflow iterator::apply( [&](auto& index){ @@ -251,7 +268,11 @@ saw::error_or step( flow_out.apply(fields,index,t_i); // equi.apply(fields,index,t_i); collision.apply(fields,macros,index,t_i); - break; + break; + case 5u: + // Corners + bb.apply(fields,index,t_i); + break; default: break; } diff --git a/lib/core/c++/boundary.hpp b/lib/core/c++/boundary.hpp index 217013d..4dbbdf8 100644 --- a/lib/core/c++/boundary.hpp +++ b/lib/core/c++/boundary.hpp @@ -174,14 +174,13 @@ class component, Encode> final { private: saw::data rho_setting_; public: - component(const saw::data& density_setting__): - rho_setting_{density_setting__} + component(const saw::data& rho_setting__): + rho_setting_{rho_setting__} {} template void apply(const saw::data& field, saw::data> index, saw::data time_step) const { using dfi = df_info; - constexpr int known_dir = Dir ? 1 : -1; bool is_even = ((time_step.get() % 2) == 0); -- cgit v1.2.3