From 530bd236a0aee9c22b87de236164c59afe6adae7 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Sat, 25 Jul 2026 23:17:18 +0200 Subject: dangling --- examples/poiseulle_particles_2d_gpu/sim.cpp | 29 ++++++++++++++++++----------- modules/core/c++/hlbm.hpp | 3 +++ 2 files changed, 21 insertions(+), 11 deletions(-) diff --git a/examples/poiseulle_particles_2d_gpu/sim.cpp b/examples/poiseulle_particles_2d_gpu/sim.cpp index 865babf..bb5787f 100644 --- a/examples/poiseulle_particles_2d_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_gpu/sim.cpp @@ -22,18 +22,11 @@ using LbmArgs = Args< } template -saw::error_or lbm_main(int argc, char** argv){ +saw::error_or lbm_main(const saw::data& args){ using namespace kel::lbm; using dfi = df_info; - auto eo_args = saw::parse_args(argc,argv); - if(eo_args.is_error()){ - return std::move(eo_args.get_error()); - } - auto& args = eo_args.get_value(); - (void)args; - auto eo_lbm_dir = output_directory(); if(eo_lbm_dir.is_error()){ return std::move(eo_lbm_dir.get_error()); @@ -198,14 +191,28 @@ saw::error_or lbm_main(int argc, char** argv){ using FloatT = kel::lbm::sch::Float32; -saw::error_or kmain(int argc, char** argv){ +saw::error_or k_main(int argc, char** argv){ using namespace kel::lbm; + + saw::data args; + { + auto eov = saw::parse_args(args,argc,argv); + return eov; + } + auto& an = args.template get<"args">(); + auto& coupling = an.template get<"coupling">(); + + if(coupling.view() == "hlbm"){ + return lbm_main(args); + //}else if(coupling.view() == "psm"){ + // return kmain(args); + } - // auto eo_args = saw::parse_args< + return saw::make_error("Invalid coupling"); } int main(int argc, char** argv){ - auto eov = lbm_main(argc, argv); + auto eov = k_main(argc, argv); if(eov.is_error()){ auto& err = eov.get_error(); std::cerr<<"[Error] "< class component final { -- cgit v1.2.3 From fb3a5f5bcba542974cf2f3f1f401ec717e662942 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Mon, 27 Jul 2026 13:32:05 +0200 Subject: Cleaning up simulation setup --- examples/poiseulle_particles_2d_gpu/common.hpp | 2 +- examples/poiseulle_particles_2d_gpu/sim.cpp | 16 ++++++++-------- examples/poiseulle_particles_2d_gpu/step.hpp | 5 +++-- 3 files changed, 12 insertions(+), 11 deletions(-) diff --git a/examples/poiseulle_particles_2d_gpu/common.hpp b/examples/poiseulle_particles_2d_gpu/common.hpp index 6c05b64..90b7687 100644 --- a/examples/poiseulle_particles_2d_gpu/common.hpp +++ b/examples/poiseulle_particles_2d_gpu/common.hpp @@ -55,7 +55,7 @@ using MacroStruct = Struct< >; template -using ParticleSpheroidGroup = ParticleGroup>; +using ParticleSpheroidGroup = ParticleGroup>; } } diff --git a/examples/poiseulle_particles_2d_gpu/sim.cpp b/examples/poiseulle_particles_2d_gpu/sim.cpp index bb5787f..a4f08f0 100644 --- a/examples/poiseulle_particles_2d_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_gpu/sim.cpp @@ -22,7 +22,7 @@ using LbmArgs = Args< } template -saw::error_or lbm_main(const saw::data& args){ +saw::error_or lbm_main(const saw::data& args){ using namespace kel::lbm; using dfi = df_info; @@ -111,7 +111,6 @@ saw::error_or lbm_main(const saw::data& args){ } } 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); @@ -188,13 +187,13 @@ saw::error_or lbm_main(const saw::data& args){ return saw::make_void(); } - using FloatT = kel::lbm::sch::Float32; +using DescT = kel::lbm::sch::D2Q9; saw::error_or k_main(int argc, char** argv){ using namespace kel::lbm; - saw::data args; + saw::data args; { auto eov = saw::parse_args(args,argc,argv); return eov; @@ -202,17 +201,18 @@ saw::error_or k_main(int argc, char** argv){ auto& an = args.template get<"args">(); auto& coupling = an.template get<"coupling">(); - if(coupling.view() == "hlbm"){ - return lbm_main(args); + if(coupling.stl_view() == "hlbm"){ + return lbm_main(args); //}else if(coupling.view() == "psm"){ - // return kmain(args); + // return kmain(args); } return saw::make_error("Invalid coupling"); } + int main(int argc, char** argv){ - auto eov = k_main(argc, argv); + auto eov = k_main(argc, argv); if(eov.is_error()){ auto& err = eov.get_error(); std::cerr<<"[Error] "< step( auto& info_f = fields.template get<"info">(); auto& porous_f = macros.template get<"porosity">(); - component> particle; + // component> particle; // auto coll_ev = q.submit([&](acpp::sycl::handler& h){ component> collision{0.8}; component> bb; - component,encode::Sycl> abb; saw::data> rho_b; rho_b.at({}) = 1.0; @@ -77,11 +76,13 @@ saw::error_or step( }); }).wait(); + /* q.submit([&](acpp::sycl::handler& h){ h.parallel_for(acpp::sycl::range<1u>{1u}, [=](acpp::sycl::id<1u> idx){ particle.apply(fields,macros,particles,{{0u}},t_i); }); }).wait(); + */ // Step /* -- cgit v1.2.3 From a826f7d9d4bb48565036dd665cb4c6a54622010d Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Mon, 27 Jul 2026 21:53:27 +0200 Subject: Dangling --- .../.nix/derivation.nix | 41 ++ .../SConscript | 34 ++ .../SConstruct | 81 ++++ .../sim.cpp | 472 +++++++++++++++++++++ examples/poiseulle_particles_2d_gpu/init.hpp | 6 +- examples/poiseulle_particles_2d_gpu/sim.cpp | 8 +- examples/poiseulle_particles_2d_gpu/step.hpp | 55 ++- modules/core/c++/fplbm.hpp | 117 +++++ 8 files changed, 779 insertions(+), 35 deletions(-) create mode 100644 examples/moving_poiseulle_particles_2d_fplbm_gpu/.nix/derivation.nix create mode 100644 examples/moving_poiseulle_particles_2d_fplbm_gpu/SConscript create mode 100644 examples/moving_poiseulle_particles_2d_fplbm_gpu/SConstruct create mode 100644 examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp diff --git a/examples/moving_poiseulle_particles_2d_fplbm_gpu/.nix/derivation.nix b/examples/moving_poiseulle_particles_2d_fplbm_gpu/.nix/derivation.nix new file mode 100644 index 0000000..672927c --- /dev/null +++ b/examples/moving_poiseulle_particles_2d_fplbm_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_fplbm_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_fplbm_gpu/SConscript b/examples/moving_poiseulle_particles_2d_fplbm_gpu/SConscript new file mode 100644 index 0000000..5e19236 --- /dev/null +++ b/examples/moving_poiseulle_particles_2d_fplbm_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_fplbm_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_fplbm_gpu/SConstruct b/examples/moving_poiseulle_particles_2d_fplbm_gpu/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/moving_poiseulle_particles_2d_fplbm_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_fplbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp new file mode 100644 index 0000000..3164945 --- /dev/null +++ b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp @@ -0,0 +1,472 @@ +#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 +using ParticleSpheroidGroup = ParticleGroup>; +} + +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}} + ); + + saw::data> rad; + rad.at({}) = dim_y * 0.1f; + saw::data> dense; + dense.at({}) = 1.0f; + particles = create_spheroid_particle_group(rad,dense,{16u}); + { + auto& rb = particles.template get<"particles">().at({0u}).template get<"rigid_body">(); + auto& rbp = rb.template get<"position">(); + rbp.at({{0u}}) = 0.25 * dim_x; + rbp.at({{1u}}) = 0.5 * dim_y; + rb.template get<"position_old">() = rbp; + saw::data> zero; + zero.at({{0u}}) = 0.0f; + zero.at({{1u}}) = 0.0f; + rb.template get<"acceleration">() = zero; + rb.template get<"rotation">() = {}; + rb.template get<"rotation_old">() = {}; + rb.template get<"angular_acceleration">() = {}; + } + + return saw::make_void(); +} + +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">(); + + q.submit([&](acpp::sycl::handler& h){ + component> hlbm_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); + }); + }).wait(); + + q.submit([&](acpp::sycl::handler& h){ + component> hlbm_one_part; + + h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){ + saw::data> index; + for(uint64_t i = 0u; i < 1u; ++i){ + index.at({{i}}).set(idx[i]); + } + + hlbm_one_part.apply(fields,macros,particles,index,t_i,{16u}); + }); + }).wait(); + + // auto coll_ev = + q.submit([&](acpp::sycl::handler& h){ + component> collision{0.8}; + component> bb; + + component,encode::Sycl> flow_in{ + [&](){ + uint64_t target_t_i = 8u; + if(t_i.get() < target_t_i){ + return 1.0 + (0.01 / target_t_i) * t_i.get(); + } + return 1.01; + }() + }; + 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); + collision.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.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>>(); + auto lbm_parts_data_ptr = saw::heap>>(); + //saw::data>,encode::Sycl>& particles + + std::cout<<"Estimated Bytes: "<,sch::MacroStruct>().get()<(*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_parts_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_parts_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; + } + } + { + auto eov = dev.copy_to_device(*lbm_parts_data_ptr,lbm_sycl_parts_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_parts_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,lsdp_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_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_parts_data,*lbm_parts_data_ptr); + if(eov.is_error()){ + return eov; + } + auto& p = lbm_parts_data_ptr->template get<"particles">().at({0u}).template get<"rigid_body">(); + auto& ppos = p.template get<"position">(); + auto& pacc = p.template get<"acceleration">(); + std::cout<<"Particle: "<> 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( +template +saw::error_or init( const converter& conv, saw::data>& fields, saw::data>& macros, @@ -117,7 +117,7 @@ saw::error_or setup_initial_conditions( radius_p.at({}).set(2); saw::data> dense_p; dense_p.at({}).set(1); - particles = create_spheroid_particle_group(radius_p, dense_p, {{16u}}); + particles = create_spheroid_particle_group(radius_p, 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 a4f08f0..fd6cdca 100644 --- a/examples/poiseulle_particles_2d_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_gpu/sim.cpp @@ -81,7 +81,7 @@ saw::error_or lbm_main(const saw::data& args){ sycl_q.wait(); { - auto eov = setup_initial_conditions(conv,*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_data_ptr); + auto eov = init(conv,*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_data_ptr); if(eov.is_error()){ return eov; } @@ -115,21 +115,21 @@ saw::error_or lbm_main(const saw::data& args){ auto lsdm_view = make_view(lbm_sycl_macro_data); auto lsdp_view = make_view(lbm_sycl_particle_data); - saw::data time_steps{16u*4096ul}; + saw::data time_steps{64u*1024ul}; auto& info_f = lsd_view.template get<"info">(); for(saw::data i{0u}; i < time_steps and krun; ++i){ // BC + Collision { - auto eov = step(conv,lsd_view,lsdm_view,lsdp_view,i,dev); + auto eov = step(conv,lsd_view,lsdm_view,lsdp_view,i,dev); if(eov.is_error()){ return eov; } } sycl_q.wait(); - if(i.get() % 32u == 0u){ + if(i.get() % 64u == 0u){ { auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr); if(eov.is_error()){ diff --git a/examples/poiseulle_particles_2d_gpu/step.hpp b/examples/poiseulle_particles_2d_gpu/step.hpp index c1b148f..37eb4e7 100644 --- a/examples/poiseulle_particles_2d_gpu/step.hpp +++ b/examples/poiseulle_particles_2d_gpu/step.hpp @@ -5,7 +5,7 @@ namespace kel { namespace lbm { -template +template saw::error_or step( const converter& conv, saw::data>,encode::Sycl>& fields, @@ -17,28 +17,44 @@ saw::error_or step( auto& q = dev.get_handle(); auto& info_f = fields.template get<"info">(); auto& porous_f = macros.template get<"porosity">(); + + q.submit([&](acpp::sycl::handler& h){ + component> hlbm_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); + }); + }).wait(); + + q.submit([&](acpp::sycl::handler& h){ + component> hlbm_one_part; + + h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){ + saw::data> index; + for(uint64_t i = 0u; i < 1u; ++i){ + index.at({{i}}).set(idx[i]); + } - // component> particle; + hlbm_one_part.apply(fields,macros,particles,index,t_i,{16u}); + }); + }).wait(); // auto coll_ev = q.submit([&](acpp::sycl::handler& h){ component> collision{0.8}; component> bb; - 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 = 16u; if(t_i.get() < target_t_i){ - return 1.0 + (0.0002 / target_t_i) * t_i.get(); + return 1.0 + (0.01 / target_t_i) * t_i.get(); } - return 1.0002; + return 1.01; }() }; component,encode::Sycl> flow_out{1.0}; @@ -62,12 +78,10 @@ saw::error_or step( 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: @@ -76,21 +90,6 @@ saw::error_or step( }); }).wait(); - /* - q.submit([&](acpp::sycl::handler& h){ - h.parallel_for(acpp::sycl::range<1u>{1u}, [=](acpp::sycl::id<1u> idx){ - particle.apply(fields,macros,particles,{{0u}},t_i); - }); - }).wait(); - */ - - // Step - /* - q.submit([&](acpp::sycl::handler& h){ - // h.depends_on(collision_ev); - }).wait(); - */ - return saw::make_void(); } diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp index c974472..34806aa 100644 --- a/modules/core/c++/fplbm.hpp +++ b/modules/core/c++/fplbm.hpp @@ -5,9 +5,31 @@ namespace kel { namespace lbm { namespace cmpt { +struct FpLbmReset{}; struct FpLbm {}; +struct FpLbmOneParticle {}; struct FpLbmOneParticleNoVelocity{}; } +namespace method { +struct FpLbm {}; +} + +template +class component final { +public: + using Component = cmpt::FpLbmReset; +private: +public: + component() = default; + + template + void apply(const saw::data& field, const saw::data& macros, saw::data> index, saw::data time_step) const { + auto& por_f = macros.template get<"porosity">(); + + auto& por = por_f.at(index); + por.at({}) = 1.0; + } +}; template class component final { @@ -84,6 +106,101 @@ public: } }; +template +class component final { +public: + using Component = cmpt::FpLbmOneParticle; +private: +public: + component() = default; + + 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">(); + auto& por_f = macros.template get<"porosity">(); + + // 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); + + auto F_i = w * force_projection; + + dfs.at({i}) = dfs.at({i}) + + frequency_ + * (eq.at(i) - dfs.at({i}) ) + + F_i.at({}) + * (saw::data{1} - saw::data{0.5f} * frequency_); + } + } + + +}; + template class component final { private: -- cgit v1.2.3