From 3077dce642f110a3e010d154a9687310ab71de43 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Sat, 20 Dec 2025 19:31:56 +0100 Subject: Working forcing out --- examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp | 8 +- examples/poiseulle_3d_gpu/.nix/derivation.nix | 39 +++ examples/poiseulle_3d_gpu/SConscript | 34 ++ examples/poiseulle_3d_gpu/SConstruct | 79 +++++ examples/poiseulle_3d_gpu/poiseulle_3d_gpu.cpp | 381 ++++++++++++++++++++ .../poiseulle_particles_2d_gpu/.nix/derivation.nix | 39 +++ examples/poiseulle_particles_2d_gpu/SConscript | 34 ++ examples/poiseulle_particles_2d_gpu/SConstruct | 79 +++++ .../poiseulle_2d_gpu.cpp | 386 +++++++++++++++++++++ 9 files changed, 1075 insertions(+), 4 deletions(-) create mode 100644 examples/poiseulle_3d_gpu/.nix/derivation.nix create mode 100644 examples/poiseulle_3d_gpu/SConscript create mode 100644 examples/poiseulle_3d_gpu/SConstruct create mode 100644 examples/poiseulle_3d_gpu/poiseulle_3d_gpu.cpp create mode 100644 examples/poiseulle_particles_2d_gpu/.nix/derivation.nix create mode 100644 examples/poiseulle_particles_2d_gpu/SConscript create mode 100644 examples/poiseulle_particles_2d_gpu/SConstruct create mode 100644 examples/poiseulle_particles_2d_gpu/poiseulle_2d_gpu.cpp diff --git a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp index 2c50dd4..89ad709 100644 --- a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp +++ b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp @@ -344,17 +344,17 @@ saw::error_or kel_main(int argc, char** argv){ } - for(uint64_t i = 0u; i < 1024u*1204u; ++i){ + uint64_t time_max = 1024u*128u; + for(uint64_t i = 0u; i < time_max; ++i){ + lbm::step(cells,macro_cells,meta,i,sycl_q); sycl_q.wait(); - if(i%128u == 0u){ + if(i%4u == 0u){ std::string vtk_f_name{"tmp/poiseulle_2d_gpu_"}; vtk_f_name += std::to_string(i) + ".vtk"; // write_vtk_file(vtk_f_name,host_cells); sycl_q.memcpy(&host_cells[0u], macro_cells, x_d * y_d * sizeof(saw::data) ).wait(); lbm::write_vtk_file(vtk_f_name, &host_cells[0], meta); } - lbm::step(cells,macro_cells,meta,i,sycl_q); - } sycl_q.wait(); diff --git a/examples/poiseulle_3d_gpu/.nix/derivation.nix b/examples/poiseulle_3d_gpu/.nix/derivation.nix new file mode 100644 index 0000000..ee78fc8 --- /dev/null +++ b/examples/poiseulle_3d_gpu/.nix/derivation.nix @@ -0,0 +1,39 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, python3 +, pname +, version +, adaptive-cpp +, kel-lbm +}: + +stdenv.mkDerivation { + pname = pname + "-examples-" + "poiseulle_3d_gpu"; + inherit version; + src = ./..; + + nativeBuildInputs = [ + scons + clang-tools + python3 + ]; + + buildInputs = [ + forstio.core + forstio.async + forstio.codec + forstio.codec-unit + forstio.remote + forstio.codec-json + adaptive-cpp + kel-lbm.core +# kel-lbm.sycl + ]; + + preferLocalBuild = true; + + outputs = [ "out" "dev" ]; +} diff --git a/examples/poiseulle_3d_gpu/SConscript b/examples/poiseulle_3d_gpu/SConscript new file mode 100644 index 0000000..2525367 --- /dev/null +++ b/examples/poiseulle_3d_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, ['poiseulle_3d_gpu.cpp'], shared=False); +examples_env.poiseulle_3d_gpu = examples_env.Program('#bin/poiseulle_3d_gpu', [examples_objects]); + +# Set Alias +env.examples = [ + examples_env.poiseulle_3d_gpu +]; +env.Alias('examples', env.examples); +env.targets += ['examples']; +env.Install('$prefix/bin/', env.examples); diff --git a/examples/poiseulle_3d_gpu/SConstruct b/examples/poiseulle_3d_gpu/SConstruct new file mode 100644 index 0000000..a7201cb --- /dev/null +++ b/examples/poiseulle_3d_gpu/SConstruct @@ -0,0 +1,79 @@ +#!/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' + ] +); +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/poiseulle_3d_gpu/poiseulle_3d_gpu.cpp b/examples/poiseulle_3d_gpu/poiseulle_3d_gpu.cpp new file mode 100644 index 0000000..c221f42 --- /dev/null +++ b/examples/poiseulle_3d_gpu/poiseulle_3d_gpu.cpp @@ -0,0 +1,381 @@ +#include +#include + +#include + +template +using SyclHostAlloc = acpp::sycl::usm_allocator, acpp::sycl::usm::alloc::host>; + +namespace kel { +namespace lbm { +namespace sch { +using namespace saw::schema; + +/** + * Basic distribution function + * Base type + * D + * Q + * Scalar factor + * D factor + * Q factor + */ +using T = Float64; +using D2Q9 = Descriptor<2u,9u>; + +using DfCell = Cell; +using CellInfo = Cell; + +/** + * Basic type for simulation + */ +using CellStruct = Struct< + Member, + Member, + Member +>; + +using MacroStruct = Struct< + Member, "velocity">, + Member +>; + +} + +namespace cmpt { +template +struct PressureBoundaryRestrictedVelocityTo {}; +} + +template +struct component> { +private: + saw::data pressure_setting_; + saw::data rho_setting_; +public: + component(const saw::data& pressure_setting__): + pressure_setting_{pressure_setting__}, + rho_setting_{pressure_setting__ * df_info::inv_cs2} + {} + + template + void apply( + saw::data* field, + saw::data> index, + const saw::data>& meta, + uint64_t time_step + ) const { + + using dfi = df_info; + + auto flat_ind = flatten_index::apply(index,meta); + + bool is_even = ((time_step % 2u) == 0u); + auto& cell = field[flat_ind.get()]; + + auto& info = cell.template get<"info">(); + if(info({0u}).get() == 0u){ + return; + } + auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + + /** + * Sum all known DFs + */ + saw::data sum_df{0u}; + for(saw::data k{0u}; k < saw::data{Desc::Q}; ++k){ + auto c_k = dfi::directions[k.get()]; + auto flat_ind_n = flatten_index::apply({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}},meta); + auto& cell_n = field[flat_ind_n.get()]; + auto& info_n = cell_n.template get<"info">(); + auto info_n_val = info_n({0u}); + auto k_opp = dfi::opposite_index[k.get()]; + + if(info_n_val.get() > 0u){ + sum_df += dfs_old({k_opp}); + } + } + /** + * Get the sum of the unknown dfs and precalculate the direction + */ + constexpr int known_dir = East ? 1 : -1; + auto sum_unknown_dfs = (rho_setting_ - sum_df) * saw::data{known_dir}; + + for(saw::data k{0u}; k < saw::data{Desc::Q}; ++k){ + auto c_k = dfi::directions[k.get()]; + auto flat_ind_n = flatten_index::apply({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}},meta); + auto& cell_n = field[flat_ind_n.get()]; + auto& info_n = cell_n.template get<"info">(); + auto info_n_val = info_n({0u}); + + if(info_n_val.get() > 0u){ + sum_unknown_dfs += dfs_old({k}) * c_k[0u]; + } + } + + auto vel_x = sum_unknown_dfs / rho_setting_; + + if constexpr (East) { + dfs_old({2u}) = dfs_old({1u}) + saw::data{2.0 / 3.0} * rho_setting_ * vel_x; + dfs_old({6u}) = dfs_old({5u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({4u}) - dfs_old({3u})); + dfs_old({8u}) = dfs_old({7u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({3u}) - dfs_old({4u})); + }else if constexpr (not East){ + dfs_old({1u}) = dfs_old({2u}) - saw::data{2.0 / 3.0} * rho_setting_ * vel_x; + dfs_old({5u}) = dfs_old({6u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({3u}) - dfs_old({4u})); + dfs_old({7u}) = dfs_old({8u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({4u}) - dfs_old({3u})); + } + } +}; + +/** + * This is massively hacky and expects a lot of conditions + * Either this or mirrored along the horizontal line works + * + * 0 - 2 - 2 + * 0 - 3 - 1 + * 0 - 3 - 1 + * ......... + * 0 - 3 - 1 + * 0 - 2 - 2 + * + */ + +template +saw::error_or set_geometry( + saw::data* cells, + const saw::data>& meta, + acpp::sycl::queue& sycl_q +){ + using namespace kel::lbm; + + saw::data rho{1.0}; + saw::data> vel{{0.0,0.0}}; + auto eq = equilibrium(rho, vel); + + sycl_q.submit([&](acpp::sycl::handler& h){ + h.parallel_for(acpp::sycl::range<2>{meta.at({0}).get(), meta.at({1}).get()},[=](acpp::sycl::id<2> idx){ + size_t i = idx[0]; + size_t j = idx[1]; + size_t acc_id = j * meta.at({0u}).get() + i; + + auto& c = cells[acc_id]; + auto& info = c.template get<"info">()({0}); + auto& dfs = c.template get<"dfs">(); + auto& dfs_old = c.template get<"dfs_old">(); + + if(i >= 2u and j >= 2u and (i+2u) < meta.at({0u}).get() and (j+2u) < meta.at({1u}).get()){ + // Fluid + info.set({2u}); + }else if(((j+2u) == meta.at({1u}).get() or j == 1u) and (i>=1u and (i+1u)= 1 and (j+1 < meta.at({1u}).get()) ) ){ + // Left input + info.set({3u}); + }else if((i+2u) == meta.at({0u}).get() and (j >= 1 and (j+1) < meta.at({1u}).get() )){ + // Right output + info.set({4u}); + }else { + info.set({0u}); + } + for(saw::data k{0u}; k < saw::data{Desc::Q}; ++k){ + dfs(k) = eq.at(k); + dfs_old(k) = eq.at(k); + } + }); + }).wait(); + + return saw::make_void(); +} + +template +void step( + saw::data* cells, + saw::data* macro_cells, + const saw::data>& meta, + uint64_t time_step, + acpp::sycl::queue& sycl_q +){ + using namespace kel::lbm; + using dfi = df_info; + + constexpr saw::data frequency{1.0 / 0.5384}; + + bool is_even = ((time_step % 2u) == 0u); + /** + * 1. Relaxation parameter \tau + */ + /* + component coll{0.5384}; + component bb; + */ + component> inlet{1.1 * dfi::cs2}; + component> outlet{1.0 * dfi::cs2}; + + + // auto collision_ev = + sycl_q.submit([&](acpp::sycl::handler& h){ + + /// Collision + h.parallel_for(acpp::sycl::range<2>{meta.at({0}).get(), meta.at({1}).get()},[=](acpp::sycl::id idx){ + size_t i = idx[0]; + size_t j = idx[1]; + size_t acc_id = j * meta.at({0u}).get() + i; + + auto& c = cells[acc_id]; + auto& info = cells[acc_id].template get<"info">(); + + switch (info({0u}).get()) { + // Bounce Back + case 1u: { + auto& dfs_old = is_even ? c.template get<"dfs_old">() : c.template get<"dfs">(); + auto df_cpy = dfs_old.copy(); + + for(uint64_t k = 1u; k < Desc::Q; ++k){ + dfs_old({k}) = df_cpy({dfi::opposite_index.at(k)}); + } + + break; + } + // Collision + case 2u: { + // coll.apply(latt_acc, {i, j}, time_step); + auto& dfs_old = is_even ? c.template get<"dfs_old">() : c.template get<"dfs">(); + + auto& macro_c = macro_cells[acc_id]; + + saw::data& rho = macro_c.template get<"pressure">(); + saw::data>& vel = macro_c.template get<"velocity">(); + + compute_rho_u(dfs_old,rho,vel); + auto eq = equilibrium(rho,vel); + + for(uint64_t k = 0u; k < Desc::Q; ++k){ + dfs_old({k}) = dfs_old({k}) + frequency * (eq.at({k}) - dfs_old({k})); + } + break; + } + case 3u: { + inlet.apply(cells, {{i,j}}, meta, time_step); + break; + } + case 4u: { + outlet.apply(cells, {{i,j}}, meta, time_step); + break; + } + default: + // Do nothing + break; + } + }); + }).wait(); + + //auto stream_ev = + sycl_q.submit([&](acpp::sycl::handler& h){ + /// Stream + h.parallel_for(acpp::sycl::range<2>{meta.at({0}).get(), meta.at({1}).get()},[=](acpp::sycl::id idx){ + size_t i = idx[0]; + size_t j = idx[1]; + size_t acc_id = j * meta.at({0u}).get() + i; + + auto& c = cells[acc_id]; + auto& info = c.template get<"info">(); + auto& dfs_new = is_even ? c.template get<"dfs">() : c.template get<"dfs_old">(); + + if (info({0u}).get() > 0u) { + for (uint64_t k = 0u; k < Desc::Q; ++k) { + auto dir = dfi::directions[dfi::opposite_index[k]]; + size_t acc_old_id = (j+dir[1]) * meta.at({0u}).get() + (i+dir[0]); + + auto& dfs_old = is_even ? cells[acc_old_id].template get<"dfs_old">() : cells[acc_old_id].template get<"dfs">(); + auto& info_old = cells[acc_old_id].template get<"info">(); + + dfs_new({k}) = dfs_old({k}); + } + } + + }); + }).wait(); +} +} +} + +template +saw::error_or kel_main(int argc, char** argv){ + using namespace kel; + + using dfi = lbm::df_info; + + auto eo_lbm_dir = lbm::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_channel_2d_gpu"; + + // Create Dir TODO + + // + lbm::converter conv { + // delta_x + {{1.0}}, + // delta_t + {{1.0}} + }; + + uint64_t x_d = 512u; + uint64_t y_d = 128u; + uint64_t z_d = 128u; + + saw::data> meta{{x_d,y_d,z_d}}; + + acpp::sycl::queue sycl_q; + SyclHostAlloc sycl_host_alloc{sycl_q}; + // SyclDeviceAlloc sycl_dev_alloc{sycl_q}; + + std::vector, SyclHostAlloc> host_cells{x_d * y_d * z_d,sycl_host_alloc}; + + saw::data* cells = acpp::sycl::malloc_device>(x_d * y_d * z_d,sycl_q); + saw::data* macro_cells = acpp::sycl::malloc_device>(x_d * y_d * z_d,sycl_q); + { + auto eov = lbm::set_geometry(cells,meta,sycl_q); + if(eov.is_error()){ + return eov; + } + } + + + for(uint64_t i = 0u; i < 1024u*128u; ++i){ + lbm::step(cells,macro_cells,meta,i,sycl_q); + sycl_q.wait(); + if(i%4u == 0u){ + std::string vtk_f_name{"tmp/poiseulle_2d_gpu_"}; + vtk_f_name += std::to_string(i) + ".vtk"; + // write_vtk_file(vtk_f_name,host_cells); + sycl_q.memcpy(&host_cells[0u], macro_cells, x_d * y_d * z_d * sizeof(saw::data) ).wait(); + lbm::write_vtk_file(vtk_f_name, &host_cells[0], meta); + } + } + + sycl_q.wait(); + acpp::sycl::free(cells, sycl_q); + acpp::sycl::free(macro_cells, sycl_q); + sycl_q.wait(); + + return saw::make_void(); +} + +int main(int argc, char** argv){ + auto eov = kel_main(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + std::cerr<<"[Error] "< 0u){ + std::cerr<<" - "< +#include + +#include + +template +using SyclHostAlloc = acpp::sycl::usm_allocator, acpp::sycl::usm::alloc::host>; + +namespace kel { +namespace lbm { +namespace sch { +using namespace saw::schema; + +/** + * Basic distribution function + * Base type + * D + * Q + * Scalar factor + * D factor + * Q factor + */ +using T = Float64; +using D2Q9 = Descriptor<2u,9u>; + +using DfCell = Cell; +using CellInfo = Cell; +using CellForceField = Cell; + +/** + * Basic type for simulation + */ +using CellStruct = Struct< + Member, + Member, + Member, + Member +>; + +using MacroStruct = Struct< + Member, "velocity">, + Member, + Member +>; + +} + +namespace cmpt { +template +struct PressureBoundaryRestrictedVelocityTo {}; +} + +template +struct component> { +private: + saw::data pressure_setting_; + saw::data rho_setting_; +public: + component(const saw::data& pressure_setting__): + pressure_setting_{pressure_setting__}, + rho_setting_{pressure_setting__ * df_info::inv_cs2} + {} + + template + void apply( + saw::data* field, + saw::data> index, + const saw::data>& meta, + uint64_t time_step + ) const { + + using dfi = df_info; + + auto flat_ind = flatten_index::apply(index,meta); + + bool is_even = ((time_step % 2u) == 0u); + auto& cell = field[flat_ind.get()]; + + auto& info = cell.template get<"info">(); + if(info({0u}).get() == 0u){ + return; + } + auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + + /** + * Sum all known DFs + */ + saw::data sum_df{0u}; + for(saw::data k{0u}; k < saw::data{Desc::Q}; ++k){ + auto c_k = dfi::directions[k.get()]; + auto flat_ind_n = flatten_index::apply({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}},meta); + auto& cell_n = field[flat_ind_n.get()]; + auto& info_n = cell_n.template get<"info">(); + auto info_n_val = info_n({0u}); + auto k_opp = dfi::opposite_index[k.get()]; + + if(info_n_val.get() > 0u){ + sum_df += dfs_old({k_opp}); + } + } + /** + * Get the sum of the unknown dfs and precalculate the direction + */ + constexpr int known_dir = East ? 1 : -1; + auto sum_unknown_dfs = (rho_setting_ - sum_df) * saw::data{known_dir}; + + for(saw::data k{0u}; k < saw::data{Desc::Q}; ++k){ + auto c_k = dfi::directions[k.get()]; + auto flat_ind_n = flatten_index::apply({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}},meta); + auto& cell_n = field[flat_ind_n.get()]; + auto& info_n = cell_n.template get<"info">(); + auto info_n_val = info_n({0u}); + + if(info_n_val.get() > 0u){ + sum_unknown_dfs += dfs_old({k}) * c_k[0u]; + } + } + + auto vel_x = sum_unknown_dfs / rho_setting_; + + if constexpr (East) { + dfs_old({2u}) = dfs_old({1u}) + saw::data{2.0 / 3.0} * rho_setting_ * vel_x; + dfs_old({6u}) = dfs_old({5u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({4u}) - dfs_old({3u})); + dfs_old({8u}) = dfs_old({7u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({3u}) - dfs_old({4u})); + }else if constexpr (not East){ + dfs_old({1u}) = dfs_old({2u}) - saw::data{2.0 / 3.0} * rho_setting_ * vel_x; + dfs_old({5u}) = dfs_old({6u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({3u}) - dfs_old({4u})); + dfs_old({7u}) = dfs_old({8u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({4u}) - dfs_old({3u})); + } + } +}; + +/** + * This is massively hacky and expects a lot of conditions + * Either this or mirrored along the horizontal line works + * + * 0 - 2 - 2 + * 0 - 3 - 1 + * 0 - 3 - 1 + * ......... + * 0 - 3 - 1 + * 0 - 2 - 2 + * + */ + +template +saw::error_or set_geometry( + saw::data* cells, + saw::data* macro_cells, + const saw::data>& meta, + acpp::sycl::queue& sycl_q +){ + using namespace kel::lbm; + + saw::data rho{1.0}; + saw::data> vel{{0.0,0.0}}; + auto eq = equilibrium(rho, vel); + + sycl_q.submit([&](acpp::sycl::handler& h){ + h.parallel_for(acpp::sycl::range<2>{meta.at({0}).get(), meta.at({1}).get()},[=](acpp::sycl::id<2> idx){ + size_t i = idx[0]; + size_t j = idx[1]; + size_t acc_id = j * meta.at({0u}).get() + i; + + auto& c = cells[acc_id]; + auto& info = c.template get<"info">()({0}); + auto& dfs = c.template get<"dfs">(); + auto& dfs_old = c.template get<"dfs_old">(); + + if(i >= 2u and j >= 2u and (i+2u) < meta.at({0u}).get() and (j+2u) < meta.at({1u}).get()){ + // Fluid + info.set({2u}); + }else if(((j+2u) == meta.at({1u}).get() or j == 1u) and (i>=1u and (i+1u)= 1 and (j+1 < meta.at({1u}).get()) ) ){ + // Left input + info.set({3u}); + }else if((i+2u) == meta.at({0u}).get() and (j >= 1 and (j+1) < meta.at({1u}).get() )){ + // Right output + info.set({4u}); + }else { + info.set({0u}); + } + for(saw::data k{0u}; k < saw::data{Desc::Q}; ++k){ + dfs(k) = eq.at(k); + dfs_old(k) = eq.at(k); + } + }); + }).wait(); + + return saw::make_void(); +} + +template +void step( + saw::data* cells, + saw::data* macro_cells, + const saw::data>& meta, + uint64_t time_step, + acpp::sycl::queue& sycl_q +){ + using namespace kel::lbm; + using dfi = df_info; + + constexpr saw::data frequency{1.0 / 0.5384}; + + bool is_even = ((time_step % 2u) == 0u); + /** + * 1. Relaxation parameter \tau + */ + /* + component coll{0.5384}; + component bb; + */ + component> inlet{1.1 * dfi::cs2}; + component> outlet{1.0 * dfi::cs2}; + + + // auto collision_ev = + sycl_q.submit([&](acpp::sycl::handler& h){ + + /// Collision + h.parallel_for(acpp::sycl::range<2>{meta.at({0}).get(), meta.at({1}).get()},[=](acpp::sycl::id idx){ + size_t i = idx[0]; + size_t j = idx[1]; + size_t acc_id = j * meta.at({0u}).get() + i; + + auto& c = cells[acc_id]; + auto& info = cells[acc_id].template get<"info">(); + + switch (info({0u}).get()) { + // Bounce Back + case 1u: { + auto& dfs_old = is_even ? c.template get<"dfs_old">() : c.template get<"dfs">(); + auto df_cpy = dfs_old.copy(); + + for(uint64_t k = 1u; k < Desc::Q; ++k){ + dfs_old({k}) = df_cpy({dfi::opposite_index.at(k)}); + } + + break; + } + // Collision + case 2u: { + // coll.apply(latt_acc, {i, j}, time_step); + auto& dfs_old = is_even ? c.template get<"dfs_old">() : c.template get<"dfs">(); + + auto& macro_c = macro_cells[acc_id]; + + saw::data& rho = macro_c.template get<"pressure">(); + saw::data>& vel = macro_c.template get<"velocity">(); + + compute_rho_u(dfs_old,rho,vel); + auto eq = equilibrium(rho,vel); + + for(uint64_t k = 0u; k < Desc::Q; ++k){ + dfs_old({k}) = dfs_old({k}) + frequency * (eq.at({k}) - dfs_old({k})); + } + break; + } + case 3u: { + inlet.apply(cells, {{i,j}}, meta, time_step); + break; + } + case 4u: { + outlet.apply(cells, {{i,j}}, meta, time_step); + break; + } + default: + // Do nothing + break; + } + }); + }).wait(); + + //auto stream_ev = + sycl_q.submit([&](acpp::sycl::handler& h){ + /// Stream + h.parallel_for(acpp::sycl::range<2>{meta.at({0}).get(), meta.at({1}).get()},[=](acpp::sycl::id idx){ + size_t i = idx[0]; + size_t j = idx[1]; + size_t acc_id = j * meta.at({0u}).get() + i; + + auto& c = cells[acc_id]; + auto& info = c.template get<"info">(); + auto& dfs_new = is_even ? c.template get<"dfs">() : c.template get<"dfs_old">(); + + if (info({0u}).get() > 1u) { + for (uint64_t k = 0u; k < Desc::Q; ++k) { + auto dir = dfi::directions[dfi::opposite_index[k]]; + size_t acc_old_id = (j+dir[1]) * meta.at({0u}).get() + (i+dir[0]); + + auto& dfs_old = is_even ? cells[acc_old_id].template get<"dfs_old">() : cells[acc_old_id].template get<"dfs">(); + auto& info_old = cells[acc_old_id].template get<"info">(); + + dfs_new({k}) = dfs_old({k}); + } + } + + }); + }).wait(); +} +} +} + +template +saw::error_or kel_main(int argc, char** argv){ + using namespace kel; + + using dfi = lbm::df_info; + + auto eo_lbm_dir = lbm::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_channel_2d_gpu"; + + // Create Dir TODO + + // + lbm::converter conv { + // delta_x + {{1.0}}, + // delta_t + {{1.0}} + }; + + uint64_t x_d = 256u; + uint64_t y_d = 64u; + + saw::data> meta{{x_d,y_d}}; + + acpp::sycl::queue sycl_q; + SyclHostAlloc sycl_host_alloc{sycl_q}; + // SyclDeviceAlloc sycl_dev_alloc{sycl_q}; + + std::vector, SyclHostAlloc> host_cells{x_d * y_d,sycl_host_alloc}; + + saw::data* cells = acpp::sycl::malloc_device>(x_d * y_d,sycl_q); + saw::data* macro_cells = acpp::sycl::malloc_device>(x_d * y_d,sycl_q); + { + auto eov = lbm::set_geometry(cells,macro_cells + ,meta,sycl_q); + if(eov.is_error()){ + return eov; + } + } + + uint64_t time_max_step = 1024u*128u; + + for(uint64_t i = 0u; i < time_max_step; ++i){ + lbm::step(cells,macro_cells,meta,i,sycl_q); + sycl_q.wait(); + if(i%1u == 4u){ + std::string vtk_f_name{"tmp/poiseulle_2d_gpu_"}; + vtk_f_name += std::to_string(i) + ".vtk"; + // write_vtk_file(vtk_f_name,host_cells); + sycl_q.memcpy(&host_cells[0u], macro_cells, x_d * y_d * sizeof(saw::data) ).wait(); + lbm::write_vtk_file(vtk_f_name, &host_cells[0], meta); + } + } + + sycl_q.wait(); + acpp::sycl::free(cells, sycl_q); + acpp::sycl::free(macro_cells, sycl_q); + sycl_q.wait(); + + return saw::make_void(); +} + +int main(int argc, char** argv){ + auto eov = kel_main(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + std::cerr<<"[Error] "< 0u){ + std::cerr<<" - "<