From 2bae8b93faf614b1acb5df7269087b4f3786da5a Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Tue, 30 Jun 2026 23:52:02 +0200 Subject: Final run --- default.nix | 2 +- examples/poiseulle_particles_2d_bgk_gpu/sim.cpp | 2 +- examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp | 49 +++++++++--------- examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp | 47 +++++++++--------- examples/poiseulle_particles_2d_psm_gpu/sim.cpp | 17 +++---- examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp | 60 ++++++++--------------- examples/stokes_drag_particle_2d_psm_gpu/sim.cpp | 39 +++++++-------- lib/core/c++/fplbm.hpp | 18 ++++++- lib/core/c++/lbm.hpp | 1 + lib/core/c++/momentum_gather.hpp | 54 ++++++++++++++++++++ lib/core/c++/particle/porosity.hpp | 8 +-- lib/core/c++/psm.hpp | 18 ++++++- 12 files changed, 185 insertions(+), 130 deletions(-) create mode 100644 lib/core/c++/momentum_gather.hpp diff --git a/default.nix b/default.nix index 02f9213..cc87f20 100644 --- a/default.nix +++ b/default.nix @@ -46,7 +46,7 @@ let src = builtins.fetchurl { url = "https://git.keldu.de/forstio-forstio/snapshot/master.tar.gz"; - sha256 = "sha256:15iqzmymza47jjx4wpc19mbg3zzwmkabpssf5y968f566n0fnb9a"; + sha256 = "sha256:0jg8a057zqqniv4cmwq0v7xmh42bn5vsjwyvpswpf9r3p5fpgxhb"; }; phases = [ "unpackPhase" "installPhase" ]; diff --git a/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp b/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp index 73964b7..ef4639e 100644 --- a/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp @@ -171,7 +171,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.65}; + component> collision{0.8}; component> bb; component,encode::Sycl> abb; diff --git a/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp b/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp index 93d4b46..3d16884 100644 --- a/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp @@ -51,6 +51,7 @@ template using MacroStruct = Struct< Member, "velocity">, Member, "density">, + Member, "force">, Member, "porosity"> >; @@ -144,6 +145,7 @@ saw::error_or setup_initial_conditions( {{1u,1u}} ); + /* iterator::apply( [&](auto& index){ saw::data> middle, ind_vec; @@ -162,6 +164,7 @@ saw::error_or setup_initial_conditions( {},// 0-index df_f.get_dims() ); + */ return saw::make_void(); } @@ -179,16 +182,20 @@ saw::error_or step( // auto coll_ev = q.submit([&](acpp::sycl::handler& h){ - component> collision{0.65}; + component> bgk{0.8}; + component> collision{0.8}; 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}; + saw::data> eps; + eps.at({}).set(1.5f); + saw::data> rad; + rad.at({}).set(dim_y*0.1); + saw::data> p_pos; + { + p_pos.at({{0u}}) = dim_x * 0.25; + p_pos.at({{1u}}) = dim_y * 0.5; + } + component> opa{p_pos,rad,eps}; component,encode::Sycl> flow_in{ [&](){ @@ -217,17 +224,18 @@ saw::error_or step( bb.apply(fields,index,t_i); break; case 2u: + opa.apply(macros,index,t_i); 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); + bgk.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); + bgk.apply(fields,macros,index,t_i); break; default: break; @@ -349,20 +357,7 @@ saw::error_or lbm_main(int argc, char** argv){ } } 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; @@ -402,6 +397,12 @@ saw::error_or lbm_main(int argc, char** argv){ return eov; } } + { + auto eov = write_csv_file(out_dir,"m",time_steps.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } sycl_q.wait(); return saw::make_void(); diff --git a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp index b09aa64..b1bd37d 100644 --- a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp @@ -48,6 +48,7 @@ template using MacroStruct = Struct< Member, "velocity">, Member, "density">, + Member, "force">, Member, "porosity"> >; @@ -140,7 +141,8 @@ saw::error_or setup_initial_conditions( df_f.get_dims(), {{1u,1u}} ); - + + /* iterator::apply( [&](auto& index){ saw::data> middle, ind_vec; @@ -159,6 +161,7 @@ saw::error_or setup_initial_conditions( {},// 0-index df_f.get_dims() ); + */ return saw::make_void(); } @@ -176,16 +179,23 @@ saw::error_or step( // auto coll_ev = q.submit([&](acpp::sycl::handler& h){ - component> collision{0.65}; + component> bgk{0.8}; + component> collision{0.8}; 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}; + saw::data> eps; + eps.at({}).set(1.5f); + saw::data> rad; + rad.at({}).set(dim_y*0.1); + saw::data> p_pos; + { + p_pos.at({{0u}}) = dim_x * 0.25; + p_pos.at({{1u}}) = dim_y * 0.5; + } + component> opa{p_pos,rad,eps}; + + component> fg; component,encode::Sycl> flow_in{ [&](){ @@ -214,32 +224,19 @@ saw::error_or step( bb.apply(fields,index,t_i); break; case 2u: + opa.apply(macros,index,t_i); 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; - } - } + fg.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); + bgk.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); + bgk.apply(fields,macros,index,t_i); break; default: break; diff --git a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp index b3ea7ae..7d3c6f5 100644 --- a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp @@ -49,6 +49,7 @@ template using MacroStruct = Struct< Member, "velocity">, Member, "density">, + Member, "force">, Member, "porosity"> >; @@ -216,17 +217,9 @@ saw::error_or step( auto& porous_f = macros.template get<"porosity">(); q.submit([&](acpp::sycl::handler& h){ + component> bgk{0.8}; component> collision{0.8}; 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}; - saw::data> eps; eps.at({}).set(1.5f); @@ -250,6 +243,7 @@ saw::error_or step( }; component,encode::Sycl> flow_out{1.0}; + component> fg; h.parallel_for(acpp::sycl::range{dim_x,dim_y}, [=](acpp::sycl::id idx){ saw::data> index; @@ -268,16 +262,17 @@ saw::error_or step( case 2u: opa.apply(macros,index,t_i); collision.apply(fields,macros,index,t_i); + fg.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); + bgk.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); + bgk.apply(fields,macros,index,t_i); break; case 5u: // Corners diff --git a/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp index a8ba21f..81498e2 100644 --- a/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp +++ b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp @@ -13,7 +13,7 @@ namespace lbm { constexpr uint64_t dim_y = 1024ul; constexpr uint64_t dim_x = dim_y * 2ul; -constexpr uint64_t particle_amount = 1ul; +// constexpr uint64_t particle_amount = 1ul; namespace sch { using namespace saw::schema; @@ -47,6 +47,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< Member, "velocity">, + Member, "force">, Member, "density">, Member, "porosity"> >; @@ -158,25 +159,6 @@ saw::error_or setup_initial_conditions( {{1u,1u}} ); - iterator::apply( - [&](auto& index){ - saw::data> middle, ind_vec; - middle.at({{0u}}) = dim_x * 0.5; - 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(); } @@ -196,12 +178,16 @@ saw::error_or step( 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}; + saw::data> eps; + eps.at({}).set(1.5f); + saw::data> rad; + rad.at({}).set(dim_y*0.1*0.25); + saw::data> p_pos; + { + p_pos.at({{0u}}) = dim_x * 0.5; + p_pos.at({{1u}}) = dim_y * 0.5; + } + component> opa{p_pos,rad,eps}; component,encode::Sycl> flow_in{ [&](){ @@ -230,6 +216,7 @@ saw::error_or step( abb.apply(fields,index,t_i); break; case 2u: + opa.apply(macros,index,t_i); collision.apply(fields,macros,index,t_i); break; case 3u: @@ -385,20 +372,7 @@ saw::error_or lbm_main(int argc, char** argv){ } } */ - 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; @@ -438,6 +412,12 @@ saw::error_or lbm_main(int argc, char** argv){ return eov; } } + { + auto eov = write_csv_file(out_dir,"m",time_steps.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } sycl_q.wait(); return saw::make_void(); diff --git a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp index 0ff6e6b..6942cc7 100644 --- a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp +++ b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp @@ -48,6 +48,7 @@ template using MacroStruct = Struct< Member, "velocity">, Member, "density">, + Member, "force">, Member, "porosity"> >; @@ -196,12 +197,16 @@ saw::error_or step( 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}; + saw::data> eps; + eps.at({}).set(1.5f); + saw::data> rad; + rad.at({}).set(dim_y*0.1*0.25); + saw::data> p_pos; + { + p_pos.at({{0u}}) = dim_x * 0.5; + p_pos.at({{1u}}) = dim_y * 0.5; + } + component> opa{p_pos,rad,eps}; component,encode::Sycl> flow_in{ [&](){ @@ -230,6 +235,7 @@ saw::error_or step( abb.apply(fields,index,t_i); break; case 2u: + opa.apply(macros,index,t_i); collision.apply(fields,macros,index,t_i); break; case 3u: @@ -385,20 +391,6 @@ saw::error_or lbm_main(int argc, char** argv){ } } */ - 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; @@ -438,7 +430,12 @@ saw::error_or lbm_main(int argc, char** argv){ return eov; } } - + { + auto eov = write_csv_file(out_dir,"m",time_steps.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } sycl_q.wait(); return saw::make_void(); } diff --git a/lib/core/c++/fplbm.hpp b/lib/core/c++/fplbm.hpp index 61def6b..c974472 100644 --- a/lib/core/c++/fplbm.hpp +++ b/lib/core/c++/fplbm.hpp @@ -114,6 +114,19 @@ public: 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; @@ -128,12 +141,13 @@ public: saw::data> min_two; min_two.at({}).set(-2); - auto force = vel * rho * min_two * por_f.at(index); + // 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){ diff --git a/lib/core/c++/lbm.hpp b/lib/core/c++/lbm.hpp index b34ec10..23d30b1 100644 --- a/lib/core/c++/lbm.hpp +++ b/lib/core/c++/lbm.hpp @@ -16,6 +16,7 @@ #include "iterator.hpp" #include "hlbm.hpp" #include "macroscopic.hpp" +#include "momentum_gather.hpp" #include "memory.hpp" #include "psm.hpp" #include "stream.hpp" diff --git a/lib/core/c++/momentum_gather.hpp b/lib/core/c++/momentum_gather.hpp new file mode 100644 index 0000000..7b1fc96 --- /dev/null +++ b/lib/core/c++/momentum_gather.hpp @@ -0,0 +1,54 @@ +#pragma once + +#include "common.hpp" +#include "macroscopic.hpp" +#include "component.hpp" +#include "equilibrium.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +// Gather the force, Luke! +struct ForceGather {}; +} + +template +class component { +private: +public: + template + void apply(const saw::data& field, const saw::data& macros, + saw::data> index, saw::data time_step) const { + + using dfi = df_info; + bool is_even = ((time_step.get() % 2) == 0); + + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + + auto& porous_f = macros.template get<"porosity">(); + auto& porous = porous_f.at(index); + + auto& force_f = macros.template get<"force">(); + auto& force = force_f.at(index); + + for(uint64_t k{0u}; k < Descriptor::D; ++k){ + force.at({{k}}).set(0); + } + + auto& dfs = dfs_old_f.at(index); + for(uint64_t i{0u}; i < Descriptor::Q; ++i){ + uint64_t i_opp = dfi::opposite_index[i]; + auto dfs_diff = dfs.at({i}) - dfs.at({i_opp}); + for(uint64_t k{0u}; k < Descriptor::D; ++k){ + force.at({{k}}) = force.at({{k}}) + dfs_diff * saw::data{dfi::directions[i][k]}; + } + } + + saw::data> one; + one.at({}) = 1.0; + auto flip_porous = one - porous; + force = force * flip_porous; + } +}; +} +} diff --git a/lib/core/c++/particle/porosity.hpp b/lib/core/c++/particle/porosity.hpp index 11185c5..2f39d2e 100644 --- a/lib/core/c++/particle/porosity.hpp +++ b/lib/core/c++/particle/porosity.hpp @@ -48,17 +48,17 @@ public: auto rad_high = rad + eps_h; if(s_dist.at({}) <= rad_low.at({})){ - por.at({}).set(1); + por.at({}).set(0); return por; } if(s_dist.at({}) >= rad_high.at({})){ - por.at({}).set(0); + por.at({}).set(1); return por; } { - typename saw::native_data_type::type inner = (std::numbers::pi / ( eps.at({}).get() )) * (s_dist - rad_low).at({}).get(); - auto cos_inner = std::cos(inner); + typename saw::native_data_type::type inner = (std::numbers::pi / ( 2 * eps.at({}).get() )) * (s_dist - rad_low).at({}).get(); + auto cos_inner = std::sin(inner); auto cos_i_2 = cos_inner * cos_inner; por.at({}).set(cos_i_2); } diff --git a/lib/core/c++/psm.hpp b/lib/core/c++/psm.hpp index 3bd3520..6dc146f 100644 --- a/lib/core/c++/psm.hpp +++ b/lib/core/c++/psm.hpp @@ -40,7 +40,7 @@ public: auto& rho_f = macros.template get<"density">(); auto& vel_f = macros.template get<"velocity">(); - + saw::data>& rho = rho_f.at(index); saw::data>& vel = vel_f.at(index); @@ -60,6 +60,22 @@ public: uint64_t i_opp = dfi::opposite_index[i]; dfs.at({i}) = dfs_cpy.at({i}) + frequency_ * (eq.at(i) - dfs_cpy.at({i})) * porous.at({}) + (dfs_cpy.at({i_opp}) - dfs_cpy.at({i}) ) * flip_porous; } + + auto& force_f = macros.template get<"force">(); + auto& force = force_f.at(index); + for(uint64_t k{0u}; k < Descriptor::D; ++k){ + force.at({{k}}).set(0); + } + + for(uint64_t i{0u}; i < Descriptor::Q; ++i){ + uint64_t i_opp = dfi::opposite_index[i]; + auto dfs_diff = dfs.at({i}) - dfs.at({i_opp}); + for(uint64_t k{0u}; k < Descriptor::D; ++k){ + force.at({{k}}) = force.at({{k}}) + dfs_diff; + } + } + + force = force * porous; } }; -- cgit v1.2.3 From 761e4912e14324ccf713b9309336816cd3283b38 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Wed, 1 Jul 2026 00:15:10 +0200 Subject: stokes drag impl and run --- default.nix | 10 + .../.nix/derivation.nix | 41 ++ .../stokes_drag_particle_2d_bgk_gpu/SConscript | 34 ++ .../stokes_drag_particle_2d_bgk_gpu/SConstruct | 81 ++++ examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp | 429 ++++++++++++++++++++ .../.nix/derivation.nix | 41 ++ .../stokes_drag_particle_2d_fplbm_gpu/SConscript | 34 ++ .../stokes_drag_particle_2d_fplbm_gpu/SConstruct | 81 ++++ examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp | 430 +++++++++++++++++++++ examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp | 58 +-- examples/stokes_drag_particle_2d_psm_gpu/sim.cpp | 60 ++- 11 files changed, 1265 insertions(+), 34 deletions(-) create mode 100644 examples/stokes_drag_particle_2d_bgk_gpu/.nix/derivation.nix create mode 100644 examples/stokes_drag_particle_2d_bgk_gpu/SConscript create mode 100644 examples/stokes_drag_particle_2d_bgk_gpu/SConstruct create mode 100644 examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp create mode 100644 examples/stokes_drag_particle_2d_fplbm_gpu/.nix/derivation.nix create mode 100644 examples/stokes_drag_particle_2d_fplbm_gpu/SConscript create mode 100644 examples/stokes_drag_particle_2d_fplbm_gpu/SConstruct create mode 100644 examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp diff --git a/default.nix b/default.nix index cc87f20..f0b4935 100644 --- a/default.nix +++ b/default.nix @@ -173,11 +173,21 @@ in rec { inherit kel; }; + stokes_drag_particle_2d_bgk_gpu = pkgs.callPackage ./examples/stokes_drag_particle_2d_bgk_gpu/.nix/derivation.nix { + inherit pname version stdenv forstio adaptive-cpp; + inherit kel; + }; + stokes_drag_particle_2d_psm_gpu = pkgs.callPackage ./examples/stokes_drag_particle_2d_psm_gpu/.nix/derivation.nix { inherit pname version stdenv forstio adaptive-cpp; inherit kel; }; + stokes_drag_particle_2d_fplbm_gpu = pkgs.callPackage ./examples/stokes_drag_particle_2d_fplbm_gpu/.nix/derivation.nix { + inherit pname version stdenv forstio adaptive-cpp; + inherit kel; + }; + stokes_drag_particle_2d_hlbm_gpu = pkgs.callPackage ./examples/stokes_drag_particle_2d_hlbm_gpu/.nix/derivation.nix { inherit pname version stdenv forstio adaptive-cpp; inherit kel; diff --git a/examples/stokes_drag_particle_2d_bgk_gpu/.nix/derivation.nix b/examples/stokes_drag_particle_2d_bgk_gpu/.nix/derivation.nix new file mode 100644 index 0000000..7ef04c4 --- /dev/null +++ b/examples/stokes_drag_particle_2d_bgk_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-" + "stokes_drag_particle_2d_bgk_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/stokes_drag_particle_2d_bgk_gpu/SConscript b/examples/stokes_drag_particle_2d_bgk_gpu/SConscript new file mode 100644 index 0000000..ff66d96 --- /dev/null +++ b/examples/stokes_drag_particle_2d_bgk_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/stokes_drag_particle_2d_bgk_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/stokes_drag_particle_2d_bgk_gpu/SConstruct b/examples/stokes_drag_particle_2d_bgk_gpu/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/stokes_drag_particle_2d_bgk_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/stokes_drag_particle_2d_bgk_gpu/sim.cpp b/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp new file mode 100644 index 0000000..911d84a --- /dev/null +++ b/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp @@ -0,0 +1,429 @@ +#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 ForceChunk = 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({}) = {0}; + 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({}) = 1.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> bgk{0.8}; + auto& collision = bgk; + component> bb; + component,encode::Sycl> abb; + + saw::data> eps; + eps.at({}).set(1.5f); + saw::data> rad; + rad.at({}).set(dim_y*0.1); + saw::data> p_pos; + { + p_pos.at({{0u}}) = dim_x * 0.25; + p_pos.at({{1u}}) = dim_y * 0.5; + } + component> opa{p_pos,rad,eps}; + + 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: + abb.apply(fields,index,t_i); + break; + case 2u: + opa.apply(macros,index,t_i); + collision.apply(fields,macros,index,t_i); + break; + case 3u: + flow_in.apply(fields,index,t_i); + // equi.apply(fields,index,t_i); + bgk.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + // equi.apply(fields,index,t_i); + bgk.apply(fields,macros,index,t_i); + break; + case 5u: + bb.apply(fields,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 / "stokes_drag_particle_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.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(); + + // 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<<" - "< +#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 ForceChunk = Chunk, 0u, dim_x, dim_y>; + +template +using MacroStruct = Struct< + Member, "velocity">, + Member, "density">, + Member, "force">, + 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({}) = {0}; + 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({}) = 1.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> bgk{0.8}; + component> collision{0.8}; + component> bb; + component,encode::Sycl> abb; + + saw::data> eps; + eps.at({}).set(1.5f); + saw::data> rad; + rad.at({}).set(dim_y*0.1); + saw::data> p_pos; + { + p_pos.at({{0u}}) = dim_x * 0.25; + p_pos.at({{1u}}) = dim_y * 0.5; + } + component> opa{p_pos,rad,eps}; + + 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: + abb.apply(fields,index,t_i); + break; + case 2u: + opa.apply(macros,index,t_i); + collision.apply(fields,macros,index,t_i); + break; + case 3u: + flow_in.apply(fields,index,t_i); + // equi.apply(fields,index,t_i); + bgk.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + // equi.apply(fields,index,t_i); + bgk.apply(fields,macros,index,t_i); + break; + case 5u: + bb.apply(fields,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 / "stokes_drag_particle_2d_fplbm_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(); + + // 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<<" - "<, 0u, dim_x, dim_y>; template using MacroStruct = Struct< Member, "velocity">, - Member, "force">, Member, "density">, + Member, "force">, Member, "porosity"> >; @@ -100,7 +100,6 @@ saw::error_or setup_initial_conditions( {{dim_x-1u,0u}}, {{dim_x, dim_y}} ); - // Overwrite with // Inflow iterator::apply( [&](auto& index){ @@ -158,7 +157,28 @@ saw::error_or setup_initial_conditions( 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(); } @@ -173,7 +193,9 @@ saw::error_or step( auto& info_f = fields.template get<"info">(); auto& porous_f = macros.template get<"porosity">(); + // auto coll_ev = q.submit([&](acpp::sycl::handler& h){ + component> bgk{0.8}; component> collision{0.8}; component> bb; component,encode::Sycl> abb; @@ -181,13 +203,15 @@ saw::error_or step( saw::data> eps; eps.at({}).set(1.5f); saw::data> rad; - rad.at({}).set(dim_y*0.1*0.25); + rad.at({}).set(dim_y*0.1/4.0); saw::data> p_pos; { p_pos.at({{0u}}) = dim_x * 0.5; p_pos.at({{1u}}) = dim_y * 0.5; } component> opa{p_pos,rad,eps}; + + component> fg; component,encode::Sycl> flow_in{ [&](){ @@ -218,21 +242,20 @@ saw::error_or step( case 2u: opa.apply(macros,index,t_i); collision.apply(fields,macros,index,t_i); + fg.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); + // equi.apply(fields,index,t_i); + bgk.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; + bgk.apply(fields,macros,index,t_i); + break; case 5u: - // Corners bb.apply(fields,index,t_i); - break; default: break; } @@ -342,9 +365,7 @@ 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); - 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){ @@ -356,8 +377,7 @@ saw::error_or lbm_main(int argc, char** argv){ } } sycl_q.wait(); - /* - if(i.get() % 32u == 0u){ + if(i.get() % 32u == 0u){ { auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr); if(eov.is_error()){ @@ -365,14 +385,12 @@ saw::error_or lbm_main(int argc, char** argv){ } } { - auto eov = write_vtk_file(out_dir,"m",i.get(), *lbm_macro_data_ptr); + 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; @@ -412,12 +430,6 @@ saw::error_or lbm_main(int argc, char** argv){ return eov; } } - { - auto eov = write_csv_file(out_dir,"m",time_steps.get(), *lbm_macro_data_ptr); - if(eov.is_error()){ - return eov; - } - } sycl_q.wait(); return saw::make_void(); diff --git a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp index 6942cc7..7a14e3f 100644 --- a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp +++ b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include @@ -13,7 +14,7 @@ namespace lbm { constexpr uint64_t dim_y = 1024ul; constexpr uint64_t dim_x = dim_y * 2ul; -constexpr uint64_t particle_amount = 1ul; +// constexpr uint64_t particle_amount = 1ul; namespace sch { using namespace saw::schema; @@ -158,11 +159,33 @@ saw::error_or setup_initial_conditions( df_f.get_dims(), {{1u,1u}} ); + + saw::data> eps; + eps.at({}).set(1.5f); + saw::data> rad; + rad.at({}).set(dim_y*0.1); + saw::data> p_pos; + { + p_pos.at({{0u}}) = dim_x * 0.25; + p_pos.at({{1u}}) = dim_y * 0.5; + } + + /* + component opa{p_pos,rad,eps}; + iterator::apply( + [&](auto& index){ + opa.apply(macros,index,{}); + }, + {},// 0-index + df_f.get_dims() + ); + */ + /* iterator::apply( [&](auto& index){ saw::data> middle, ind_vec; - middle.at({{0u}}) = dim_x * 0.5; + middle.at({{0u}}) = dim_x * 0.25; middle.at({{1u}}) = dim_y * 0.5; ind_vec.at({{0u}}) = index.at({{0u}}).template cast_to(); @@ -177,6 +200,7 @@ saw::error_or setup_initial_conditions( {},// 0-index df_f.get_dims() ); + */ return saw::make_void(); } @@ -193,6 +217,7 @@ saw::error_or step( auto& porous_f = macros.template get<"porosity">(); q.submit([&](acpp::sycl::handler& h){ + component> bgk{0.8}; component> collision{0.8}; component> bb; component,encode::Sycl> abb; @@ -200,7 +225,7 @@ saw::error_or step( saw::data> eps; eps.at({}).set(1.5f); saw::data> rad; - rad.at({}).set(dim_y*0.1*0.25); + rad.at({}).set(dim_y*0.1 / 4.0); saw::data> p_pos; { p_pos.at({{0u}}) = dim_x * 0.5; @@ -219,6 +244,7 @@ saw::error_or step( }; component,encode::Sycl> flow_out{1.0}; + component> fg; h.parallel_for(acpp::sycl::range{dim_x,dim_y}, [=](acpp::sycl::id idx){ saw::data> index; @@ -237,16 +263,17 @@ saw::error_or step( case 2u: opa.apply(macros,index,t_i); collision.apply(fields,macros,index,t_i); + fg.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); + bgk.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); + bgk.apply(fields,macros,index,t_i); break; case 5u: // Corners @@ -391,6 +418,22 @@ saw::error_or lbm_main(int argc, char** argv){ } } */ + /* + 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; @@ -430,12 +473,7 @@ saw::error_or lbm_main(int argc, char** argv){ return eov; } } - { - auto eov = write_csv_file(out_dir,"m",time_steps.get(), *lbm_macro_data_ptr); - if(eov.is_error()){ - return eov; - } - } + sycl_q.wait(); return saw::make_void(); } -- cgit v1.2.3