From bff15d076515b93621a66fe93f58d4d3221c23c7 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Wed, 12 Aug 2026 09:49:29 +0200 Subject: Preparing cleanup of boundarz and yesterdays changes --- modules/core/c++/boundary.hpp | 58 ++++++++++++++++++++++++++++++++++-- modules/core/c++/boundary/common.hpp | 8 +++++ modules/core/c++/fplbm.hpp | 2 +- 3 files changed, 64 insertions(+), 4 deletions(-) create mode 100644 modules/core/c++/boundary/common.hpp diff --git a/modules/core/c++/boundary.hpp b/modules/core/c++/boundary.hpp index 22a6f34..ae53c04 100644 --- a/modules/core/c++/boundary.hpp +++ b/modules/core/c++/boundary.hpp @@ -169,8 +169,8 @@ public: * 0 - 2 - 2 * */ -template -class component, Encode> final { +template +class component, cmpt::ZouHeHorizontal, Encode> final { private: saw::data rho_setting_; public: @@ -179,7 +179,7 @@ public: {} template - void apply(const saw::data& field, saw::data> index, saw::data time_step) const { + void apply(const saw::data& field, saw::data> index, saw::data time_step) const { using dfi = df_info; bool is_even = ((time_step.get() % 2) == 0); @@ -207,6 +207,58 @@ public: return {}; }(); + // static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly"); + + if constexpr (Dir) { + dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data{2.0 / 3.0} * rho_vel_x; + dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data{1.0 / 6.0} * rho_vel_x + saw::data{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); + dfs_old.at({8u}) = dfs_old.at({7u}) + saw::data{1.0 / 6.0} * rho_vel_x + saw::data{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); + }else if constexpr (not Dir){ + dfs_old.at({1u}) = dfs_old.at({2u}) - saw::data{2.0 / 3.0} * rho_vel_x; + dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data{1.0 / 6.0} * rho_vel_x + saw::data{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); + dfs_old.at({7u}) = dfs_old.at({8u}) - saw::data{1.0 / 6.0} * rho_vel_x + saw::data{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); + } + } +}; + +template +class component, cmpt::ZouHeHorizontal, Encode> final { +private: + saw::data rho_setting_; +public: + component(const saw::data& rho_setting__): + rho_setting_{rho_setting__} + {} + + template + void apply(const saw::data& field, saw::data> index, saw::data time_step) const { + using dfi = df_info; + + bool is_even = ((time_step.get() % 2) == 0); + + auto& info_f = field.template get<"info">(); + + // auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">(); + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + + auto& dfs_old = dfs_old_f.at(index); + + auto rho_vel_x = [&]() -> saw::data { + if constexpr (Dir){ + auto S = dfs_old.at({0u}) + + dfs_old.at({3u}) + dfs_old.at({4u}) + + (dfs_old.at({1u}) + dfs_old.at({5u}) + dfs_old.at({7u})) * 2u; + return rho_setting_ - S; + } + else if constexpr (not Dir) { + auto S = dfs_old.at({0u}) + + dfs_old.at({3u}) + dfs_old.at({4u}) + + (dfs_old.at({2u}) + dfs_old.at({6u}) + dfs_old.at({8u})) * 2u; + return S - rho_setting_; + } + return {}; + }(); + static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly"); if constexpr (Dir) { diff --git a/modules/core/c++/boundary/common.hpp b/modules/core/c++/boundary/common.hpp new file mode 100644 index 0000000..3feb451 --- /dev/null +++ b/modules/core/c++/boundary/common.hpp @@ -0,0 +1,8 @@ +#pragma once + +#include "../common.hpp" + +namespace kel { +namespace lbm { +} +} diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp index 37074b0..e2527ce 100644 --- a/modules/core/c++/fplbm.hpp +++ b/modules/core/c++/fplbm.hpp @@ -204,7 +204,7 @@ public: // vel_s is technically time the density of the particle? - force = ( vel_s * rho - vel * rho ) * two * flip_por; + force = ( vel_s * rho - vel * rho ) /* two */ * flip_por; force_p = force_p - force; }, aabb.template get<"a">(), aabb.template get<"b">()); -- cgit v1.2.3 From 919f5a625c5efeea94c0dd0de5039a27e77d3fac Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Wed, 12 Aug 2026 13:33:42 +0200 Subject: Fixing Pressure ZouHe, but also forgot that I wanted to move to momentum --- modules/core/c++/boundary.hpp | 505 ++++++++++++++++++++++++-- modules/core/c++/boundary/zou_he_pressure.hpp | 9 + 2 files changed, 486 insertions(+), 28 deletions(-) create mode 100644 modules/core/c++/boundary/zou_he_pressure.hpp diff --git a/modules/core/c++/boundary.hpp b/modules/core/c++/boundary.hpp index ae53c04..7f185ad 100644 --- a/modules/core/c++/boundary.hpp +++ b/modules/core/c++/boundary.hpp @@ -171,6 +171,8 @@ public: */ template class component, cmpt::ZouHeHorizontal, Encode> final { +public: + using Descriptor = sch::Descriptor<2u,9u>; private: saw::data rho_setting_; public: @@ -221,8 +223,10 @@ public: } }; -template -class component, cmpt::ZouHeHorizontal, Encode> final { +template +class component, cmpt::ZouHeHorizontal, Encode> final { +public: + using Descriptor = sch::Descriptor<3u,27u>; private: saw::data rho_setting_; public: @@ -240,37 +244,482 @@ public: // auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">(); auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + auto& f = dfs_old_f.at(index); + + /* + * D3Q27 indexing + * + * x=-1 x=0 x=+1 + * + * z=+1: 19 22 25 18 21 24 20 23 26 + * + * z= 0: 1 4 7 0 3 6 2 5 8 + * + * z=-1: 10 13 16 9 12 15 11 14 17 + * + * Boundary normal is x. + * + * East: + * known: cx = +1 + * unknown: cx = -1 + * + * West: + * known: cx = -1 + * unknown: cx = +1 + */ + + + /* + * ============================================================ + * rho * ux + * ============================================================ + * + * rho = S + 2 * sum(unknown) + * + * Therefore: + * + * East: + * rho*ux = rho - S + * + * West: + * rho*ux = S - rho + */ + + saw::data S; + + if constexpr (East) { + + S = + // cx = 0 + f.at({0u}) + + f.at({3u}) + + f.at({6u}) + + f.at({9u}) + + f.at({12u}) + + f.at({15u}) + + f.at({18u}) + + f.at({21u}) + + f.at({24u}) + + // cx = +1 + + saw::data{2.0} * ( + f.at({2u}) + + f.at({5u}) + + f.at({8u}) + + f.at({11u}) + + f.at({14u}) + + f.at({17u}) + + f.at({20u}) + + f.at({23u}) + + f.at({26u}) + ); + + } else { + + S = + // cx = 0 + f.at({0u}) + + f.at({3u}) + + f.at({6u}) + + f.at({9u}) + + f.at({12u}) + + f.at({15u}) + + f.at({18u}) + + f.at({21u}) + + f.at({24u}) + + // cx = -1 + + saw::data{2.0} * ( + f.at({1u}) + + f.at({4u}) + + f.at({7u}) + + f.at({10u}) + + f.at({13u}) + + f.at({16u}) + + f.at({19u}) + + f.at({22u}) + + f.at({25u}) + ); + } + + saw::data rho_ux; + + if constexpr (East) { + rho_ux = rho_setting_ - S; + } else { + rho_ux = S - rho_setting_; + } + + + /* + * ============================================================ + * Tangential momentum: rho * uy + * ============================================================ + * + * For a pressure boundary, uy and uz are obtained from the + * known populations. + * + * rho*uy: + * + * East: + * + * x=0 contribution + * + 2*x=+1 contribution + * + * West: + * + * x=0 contribution + * + 2*x=-1 contribution + */ + + saw::data rho_uy; + saw::data rho_uz; + + if constexpr (East) { + + rho_uy = + // cx = 0 + -f.at({3u}) + + f.at({6u}) + -f.at({12u}) + + f.at({15u}) + -f.at({21u}) + + f.at({24u}) + + // cx = +1 + + saw::data{2.0} * ( + -f.at({5u}) + + f.at({8u}) + -f.at({14u}) + + f.at({17u}) + -f.at({23u}) + + f.at({26u}) + ); + + rho_uz = + // cx = 0 + -f.at({9u}) + -f.at({12u}) + -f.at({15u}) + + f.at({18u}) + + f.at({21u}) + + f.at({24u}) + + // cx = +1 + + saw::data{2.0} * ( + -f.at({11u}) + -f.at({14u}) + -f.at({17u}) + +f.at({20u}) + +f.at({23u}) + +f.at({26u}) + ); + + } else { + + rho_uy = + // cx = 0 + -f.at({3u}) + + f.at({6u}) + - f.at({12u}) + + f.at({15u}) + - f.at({21u}) + + f.at({24u}) + + // cx = -1 + + saw::data{2.0} * ( + -f.at({4u}) + + f.at({7u}) + -f.at({13u}) + + f.at({16u}) + -f.at({22u}) + + f.at({25u}) + ); + + rho_uz = + // cx = 0 + -f.at({9u}) + -f.at({12u}) + -f.at({15u}) + +f.at({18u}) + +f.at({21u}) + +f.at({24u}) + + // cx = -1 + + saw::data{2.0} * ( + -f.at({10u}) + -f.at({13u}) + -f.at({16u}) + +f.at({19u}) + +f.at({22u}) + +f.at({25u}) + ); + } + + + /* + * ============================================================ + * Velocity + * ============================================================ + */ + + saw::data ux = rho_ux / rho_setting_; + saw::data uy = rho_uy / rho_setting_; + saw::data uz = rho_uz / rho_setting_; + + + /* + * ============================================================ + * Reconstruct unknown distributions + * ============================================================ + * + * Zou/He non-equilibrium bounce-back: + * + * f_i = f_opp + * + (f_eq_i - f_eq_opp) + * + * For D3Q27: + * + * axis: + * + * 6*w*rho*u = 1/6 * rho*u + * + * face diagonal: + * + * 6*w*rho*(c.u) + * = 1/9 * rho*(c.u) + * + * body diagonal: + * + * 6*w*rho*(c.u) + * = 1/36 * rho*(c.u) + * + * The factors below therefore use: + * + * 1/6 + * 1/9 + * 1/36 + * + * ============================================================ + */ + + if constexpr (East) { + + /* + * -------------------------------------------------------- + * (-1, 0, 0) <- (+1, 0, 0) + * -------------------------------------------------------- + */ + + f.at({1u}) = + f.at({2u}) + - saw::data{1.0 / 6.0} * rho_ux; + + + /* + * -------------------------------------------------------- + * (-1, -1, 0) <- (+1, -1, 0) + * -------------------------------------------------------- + */ + + f.at({4u}) = + f.at({5u}) + - saw::data{1.0 / 9.0} + * (rho_ux + rho_uy); + + + /* + * -------------------------------------------------------- + * (-1, +1, 0) <- (+1, +1, 0) + * -------------------------------------------------------- + */ + + f.at({7u}) = + f.at({8u}) + + saw::data{1.0 / 9.0} + * (-rho_ux + rho_uy); + + + /* + * -------------------------------------------------------- + * (-1, 0, -1) <- (+1, 0, -1) + * -------------------------------------------------------- + */ + + f.at({10u}) = + f.at({11u}) + - saw::data{1.0 / 9.0} + * (rho_ux + rho_uz); + + + /* + * -------------------------------------------------------- + * (-1, -1, -1) <- (+1, -1, -1) + * -------------------------------------------------------- + */ + + f.at({13u}) = + f.at({14u}) + - saw::data{1.0 / 36.0} + * (rho_ux + rho_uy + rho_uz); + + + /* + * -------------------------------------------------------- + * (-1, +1, -1) <- (+1, +1, -1) + * -------------------------------------------------------- + */ + + f.at({16u}) = + f.at({17u}) + + saw::data{1.0 / 36.0} + * (-rho_ux + rho_uy + rho_uz); + + + /* + * -------------------------------------------------------- + * (-1, 0, +1) <- (+1, 0, +1) + * -------------------------------------------------------- + */ + + f.at({19u}) = + f.at({20u}) + + saw::data{1.0 / 9.0} + * (-rho_ux + rho_uz); + + + /* + * -------------------------------------------------------- + * (-1, -1, +1) <- (+1, -1, +1) + * -------------------------------------------------------- + */ + + f.at({22u}) = + f.at({23u}) + + saw::data{1.0 / 36.0} + * (-rho_ux - rho_uy + rho_uz); - auto& dfs_old = dfs_old_f.at(index); - auto rho_vel_x = [&]() -> saw::data { - if constexpr (Dir){ - auto S = dfs_old.at({0u}) - + dfs_old.at({3u}) + dfs_old.at({4u}) - + (dfs_old.at({1u}) + dfs_old.at({5u}) + dfs_old.at({7u})) * 2u; - return rho_setting_ - S; - } - else if constexpr (not Dir) { - auto S = dfs_old.at({0u}) - + dfs_old.at({3u}) + dfs_old.at({4u}) - + (dfs_old.at({2u}) + dfs_old.at({6u}) + dfs_old.at({8u})) * 2u; - return S - rho_setting_; - } - return {}; - }(); + /* + * -------------------------------------------------------- + * (-1, +1, +1) <- (+1, +1, +1) + * -------------------------------------------------------- + */ + + f.at({25u}) = + f.at({26u}) + + saw::data{1.0 / 36.0} + * (-rho_ux + rho_uy + rho_uz); - static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly"); + } else { - if constexpr (Dir) { - dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data{2.0 / 3.0} * rho_vel_x; - dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data{1.0 / 6.0} * rho_vel_x + saw::data{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); - dfs_old.at({8u}) = dfs_old.at({7u}) + saw::data{1.0 / 6.0} * rho_vel_x + saw::data{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); - }else if constexpr (not Dir){ - dfs_old.at({1u}) = dfs_old.at({2u}) - saw::data{2.0 / 3.0} * rho_vel_x; - dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data{1.0 / 6.0} * rho_vel_x + saw::data{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); - dfs_old.at({7u}) = dfs_old.at({8u}) - saw::data{1.0 / 6.0} * rho_vel_x + saw::data{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); + /* + * -------------------------------------------------------- + * (+1, 0, 0) <- (-1, 0, 0) + * -------------------------------------------------------- + */ + + f.at({2u}) = + f.at({1u}) + + saw::data{1.0 / 6.0} * rho_ux; + + + /* + * -------------------------------------------------------- + * (+1, -1, 0) <- (-1, -1, 0) + * -------------------------------------------------------- + */ + + f.at({5u}) = + f.at({4u}) + + saw::data{1.0 / 9.0} + * (rho_ux - rho_uy); + + + /* + * -------------------------------------------------------- + * (+1, +1, 0) <- (-1, +1, 0) + * -------------------------------------------------------- + */ + + f.at({8u}) = + f.at({7u}) + + saw::data{1.0 / 9.0} + * (rho_ux + rho_uy); + + + /* + * -------------------------------------------------------- + * (+1, 0, -1) <- (-1, 0, -1) + * -------------------------------------------------------- + */ + + f.at({11u}) = + f.at({10u}) + + saw::data{1.0 / 9.0} + * (rho_ux - rho_uz); + + + /* + * -------------------------------------------------------- + * (+1, -1, -1) <- (-1, -1, -1) + * -------------------------------------------------------- + */ + + f.at({14u}) = + f.at({13u}) + + saw::data{1.0 / 36.0} + * (rho_ux - rho_uy - rho_uz); + + + /* + * -------------------------------------------------------- + * (+1, +1, -1) <- (-1, +1, -1) + * -------------------------------------------------------- + */ + + f.at({17u}) = + f.at({16u}) + + saw::data{1.0 / 36.0} + * (rho_ux + rho_uy - rho_uz); + + + /* + * -------------------------------------------------------- + * (+1, 0, +1) <- (-1, 0, +1) + * -------------------------------------------------------- + */ + + f.at({20u}) = + f.at({19u}) + + saw::data{1.0 / 9.0} + * (rho_ux + rho_uz); + + + /* + * -------------------------------------------------------- + * (+1, -1, +1) <- (-1, -1, +1) + * -------------------------------------------------------- + */ + + f.at({23u}) = + f.at({22u}) + + saw::data{1.0 / 36.0} + * (rho_ux - rho_uy + rho_uz); + + + /* + * -------------------------------------------------------- + * (+1, +1, +1) <- (-1, +1, +1) + * -------------------------------------------------------- + */ + + f.at({26u}) = + f.at({25u}) + + saw::data{1.0 / 36.0} + * (rho_ux + rho_uy + rho_uz); + } } - } }; diff --git a/modules/core/c++/boundary/zou_he_pressure.hpp b/modules/core/c++/boundary/zou_he_pressure.hpp new file mode 100644 index 0000000..7fc8591 --- /dev/null +++ b/modules/core/c++/boundary/zou_he_pressure.hpp @@ -0,0 +1,9 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +} +} -- cgit v1.2.3 From 7220084ddae1f4c6007f5111913cb7f87a794657 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Thu, 13 Aug 2026 13:11:05 +0200 Subject: Renaming velocity to momentum. for now --- examples/cavity_2d/cavity_2d.cpp | 4 +- examples/cavity_2d_gpu/cavity_2d_gpu.cpp | 4 +- examples/cavity_3d.cpp | 4 +- examples/meta/meta.cpp | 6 +- .../sim.cpp | 4 +- .../moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp | 4 +- .../moving_poiseulle_particles_2d_psm_gpu/sim.cpp | 4 +- examples/poiseulle_2d/poiseulle_2d.cpp | 4 +- examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp | 8 +- .../poiseulle_moving_particle_2d_psm_gpu/sim.cpp | 4 +- examples/poiseulle_particles_2d_bgk_gpu/sim.cpp | 4 +- examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp | 4 +- examples/poiseulle_particles_2d_gpu/common.hpp | 2 +- examples/poiseulle_particles_2d_gpu/init.hpp | 2 +- examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp | 4 +- examples/poiseulle_particles_2d_ibm_gpu/sim.cpp | 6 +- examples/poiseulle_particles_2d_psm_gpu/sim.cpp | 4 +- .../poiseulle_particles_channel_2d.cpp | 10 +- .../.nix/derivation.nix | 44 ++ .../SConscript | 34 ++ .../SConstruct | 81 ++++ .../common.hpp | 67 +++ .../init.hpp | 127 +++++ .../schaefer_turek_durst_krause_rannbacher/sim.cpp | 226 +++++++++ .../step.hpp | 97 ++++ examples/settling_cubes_2d_ibm_gpu/sim.cpp | 8 +- examples/settling_spheres_2d_hlbm_gpu/sim.cpp | 6 +- examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp | 4 +- examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp | 4 +- examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp | 4 +- examples/stokes_drag_particle_2d_psm_gpu/sim.cpp | 4 +- modules/core/c++/boundary.hpp | 36 +- modules/core/c++/collision.hpp | 4 +- modules/core/c++/component/psm.hpp | 2 +- modules/core/c++/fplbm.hpp | 8 +- modules/core/c++/hlbm.hpp | 4 +- modules/core/c++/macroscopic.hpp | 1 + modules/core/c++/particle.hpp | 2 +- modules/core/c++/psm.hpp | 516 ++++++++++++++++----- modules/core/c++/write_csv.hpp | 2 +- modules/core/c++/write_vtk.hpp | 2 +- modules/core/tests/equilibrium.cpp | 2 +- modules/core/tests/vtk_write.cpp | 10 +- poc/triangle_stream/common.hpp | 2 +- 44 files changed, 1179 insertions(+), 200 deletions(-) create mode 100644 examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix create mode 100644 examples/schaefer_turek_durst_krause_rannbacher/SConscript create mode 100644 examples/schaefer_turek_durst_krause_rannbacher/SConstruct create mode 100644 examples/schaefer_turek_durst_krause_rannbacher/common.hpp create mode 100644 examples/schaefer_turek_durst_krause_rannbacher/init.hpp create mode 100644 examples/schaefer_turek_durst_krause_rannbacher/sim.cpp create mode 100644 examples/schaefer_turek_durst_krause_rannbacher/step.hpp diff --git a/examples/cavity_2d/cavity_2d.cpp b/examples/cavity_2d/cavity_2d.cpp index 702661a..af22796 100644 --- a/examples/cavity_2d/cavity_2d.cpp +++ b/examples/cavity_2d/cavity_2d.cpp @@ -42,7 +42,7 @@ using CellStruct = Struct< template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member >; @@ -406,7 +406,7 @@ int main(){ auto dfs = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); auto& rho = macros.at({{i,j}}).template get<"pressure">(); - auto& vel = macros.at({{i,j}}).template get<"velocity">(); + auto& vel = macros.at({{i,j}}).template get<"momentum">(); compute_rho_u(dfs,rho,vel); }, lattice); diff --git a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp index 0fc490a..fe44baf 100644 --- a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp +++ b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp @@ -47,7 +47,7 @@ using CellStruct = CellFieldStruct< template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member >; @@ -362,7 +362,7 @@ int main(){ auto dfs_field = even_step ? dfs_old : dfs; auto& rho = macros.at(index).template get<"pressure">(); - auto& vel = macros.at(index).template get<"velocity">(); + auto& vel = macros.at(index).template get<"momentum">(); compute_rho_u(dfs_field[j],rho,vel); }, {{0u,0u}}, meta); std::string vtk_f_name{"tmp/cavity_2d_gpu_"}; diff --git a/examples/cavity_3d.cpp b/examples/cavity_3d.cpp index 64356be..5ae4834 100644 --- a/examples/cavity_3d.cpp +++ b/examples/cavity_3d.cpp @@ -46,7 +46,7 @@ using CellStruct = Struct< template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member >; @@ -411,7 +411,7 @@ int main(){ auto dfs = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); auto& rho = macros.at({{i,j}}).template get<"pressure">(); - auto& vel = macros.at({{i,j}}).template get<"velocity">(); + auto& vel = macros.at({{i,j}}).template get<"momentum">(); compute_rho_u(dfs,rho,vel); }, lattice); diff --git a/examples/meta/meta.cpp b/examples/meta/meta.cpp index a5e2ae5..7232c61 100644 --- a/examples/meta/meta.cpp +++ b/examples/meta/meta.cpp @@ -75,7 +75,7 @@ saw::error_or meta_main(int argc, char** argv){ meta_lbm_info>(args_str); break; default: - return saw::make_error("Second positional needs to be a valid velocity set"); + return saw::make_error("Second positional needs to be a valid momentum set"); } } break; @@ -89,7 +89,7 @@ saw::error_or meta_main(int argc, char** argv){ meta_lbm_info>(args_str); break; default: - return saw::make_error("Second positional needs to be a valid velocity set"); + return saw::make_error("Second positional needs to be a valid momentum set"); } } break; @@ -108,7 +108,7 @@ saw::error_or meta_main(int argc, char** argv){ meta_lbm_info>(args_str); break; default: - return saw::make_error("Second positional needs to be a valid velocity set"); + return saw::make_error("Second positional needs to be a valid momentum set"); } } break; diff --git a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp index cc7649a..4c9c8bd 100644 --- a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp +++ b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp @@ -55,7 +55,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "porosity">, Member, "force"> @@ -119,7 +119,7 @@ saw::error_or init( auto& dfs_old_f = fields.template get<"dfs_old">(); auto& dfs_f = fields.template get<"dfs">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); auto& force_f = macros.template get<"force">(); diff --git a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp index 9e812e2..bd754bf 100644 --- a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp +++ b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp @@ -46,7 +46,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "porosity"> >; @@ -109,7 +109,7 @@ saw::error_or setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator::apply( diff --git a/examples/moving_poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_psm_gpu/sim.cpp index 2c4c9d1..2daa452 100644 --- a/examples/moving_poiseulle_particles_2d_psm_gpu/sim.cpp +++ b/examples/moving_poiseulle_particles_2d_psm_gpu/sim.cpp @@ -46,7 +46,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "porosity">, Member, "force"> @@ -110,7 +110,7 @@ saw::error_or init( auto& dfs_old_f = fields.template get<"dfs_old">(); auto& dfs_f = fields.template get<"dfs">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); auto& force_f = macros.template get<"force">(); diff --git a/examples/poiseulle_2d/poiseulle_2d.cpp b/examples/poiseulle_2d/poiseulle_2d.cpp index 3bc7a6c..f9db3df 100644 --- a/examples/poiseulle_2d/poiseulle_2d.cpp +++ b/examples/poiseulle_2d/poiseulle_2d.cpp @@ -41,7 +41,7 @@ using CellStruct = Struct< template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member >; @@ -375,7 +375,7 @@ int main(int argc, char** argv){ auto& dfs = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">(); auto& rho = macros.at(index).template get<"pressure">(); - auto& vel = macros.at(index).template get<"velocity">(); + auto& vel = macros.at(index).template get<"momentum">(); compute_rho_u(dfs,rho,vel); rho = rho * saw::data{dfi::cs2}; diff --git a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp index 1aaa5c0..9963c59 100644 --- a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp +++ b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp @@ -36,7 +36,7 @@ using CellStruct = Struct< >; using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member >; @@ -245,7 +245,7 @@ void step( auto& macro_c = macro_cells[acc_id]; saw::data& rho = macro_c.template get<"pressure">(); - saw::data>& vel = macro_c.template get<"velocity">(); + saw::data>& vel = macro_c.template get<"momentum">(); compute_rho_u(dfs_old,rho,vel); auto eq = equilibrium(rho,vel); @@ -262,7 +262,7 @@ void step( auto& macro_c = macro_cells[acc_id]; saw::data& rho = macro_c.template get<"pressure">(); - saw::data>& vel = macro_c.template get<"velocity">(); + saw::data>& vel = macro_c.template get<"momentum">(); compute_rho_u(dfs_old,rho,vel); auto eq = equilibrium(rho,vel); @@ -279,7 +279,7 @@ void step( auto& macro_c = macro_cells[acc_id]; saw::data& rho = macro_c.template get<"pressure">(); - saw::data>& vel = macro_c.template get<"velocity">(); + saw::data>& vel = macro_c.template get<"momentum">(); compute_rho_u(dfs_old,rho,vel); auto eq = equilibrium(rho,vel); diff --git a/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp b/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp index 1693c5e..1af0023 100644 --- a/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp +++ b/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp @@ -47,7 +47,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "porosity"> >; @@ -123,7 +123,7 @@ saw::error_or setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator::apply( diff --git a/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp b/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp index 84ce22b..9ed1646 100644 --- a/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp @@ -45,7 +45,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density"> >; @@ -102,7 +102,7 @@ saw::error_or setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); iterator::apply( [&](auto& index){ diff --git a/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp b/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp index c6aac92..f45e10e 100644 --- a/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp @@ -49,7 +49,7 @@ using ForceChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "force">, Member, "porosity"> @@ -109,7 +109,7 @@ saw::error_or setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator::apply( diff --git a/examples/poiseulle_particles_2d_gpu/common.hpp b/examples/poiseulle_particles_2d_gpu/common.hpp index 90b7687..b321344 100644 --- a/examples/poiseulle_particles_2d_gpu/common.hpp +++ b/examples/poiseulle_particles_2d_gpu/common.hpp @@ -49,7 +49,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "porosity"> >; diff --git a/examples/poiseulle_particles_2d_gpu/init.hpp b/examples/poiseulle_particles_2d_gpu/init.hpp index cdce25e..b7fdeb0 100644 --- a/examples/poiseulle_particles_2d_gpu/init.hpp +++ b/examples/poiseulle_particles_2d_gpu/init.hpp @@ -57,7 +57,7 @@ saw::error_or init( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator::apply( diff --git a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp index 6cadbd0..0061fa2 100644 --- a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp @@ -46,7 +46,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "force">, Member, "porosity"> @@ -106,7 +106,7 @@ saw::error_or setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator::apply( diff --git a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp index 640c9ab..d63515d 100644 --- a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp @@ -46,7 +46,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "force"> >; @@ -113,7 +113,7 @@ saw::error_or setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& force_f = macros.template get<"force">(); iterator::apply( @@ -218,7 +218,7 @@ saw::error_or step( q.submit([&](acpp::sycl::handler& h){ h.parallel_for(acpp::sycl::range<1u>{1u}, [=](acpp::sycl::id<1u> idx){ - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& dense_f = macros.template get<"density">(); auto& ps = particles; diff --git a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp index c2ab675..760f70a 100644 --- a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp @@ -47,7 +47,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "force">, Member, "porosity"> @@ -124,7 +124,7 @@ saw::error_or setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator::apply( diff --git a/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp b/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp index 1e6e75f..24832f2 100644 --- a/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp +++ b/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp @@ -68,7 +68,7 @@ using CellStruct = Struct< template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, Member, Member, "force"> @@ -369,7 +369,7 @@ void couple_particles_to_lattice( // Interpolate this from close U cells. // For now pick the closest U - auto& closest_u = macros.at(p_cell_pos).template get<"velocity">(); + auto& closest_u = macros.at(p_cell_pos).template get<"momentum">(); // auto p_shift = closest_u; // p_pos - p_pos_old @@ -380,7 +380,7 @@ void couple_particles_to_lattice( p_vel_rel.at({{1u}}) = p_vel_rel.at({{1u}}) / p_mass; } - // Add relative velocity to particle acceleration (Since our timestep is 1 we don't need to do anything else) + // Add relative momentum to particle acceleration (Since our timestep is 1 we don't need to do anything else) p_acc = p_acc + p_vel_rel; @@ -463,7 +463,7 @@ void couple_particles_to_lattice( } - /* So I want to include the relative velocity from both particles, but this is a bit hard considering I just assumed 0 velocity from the other party + /* So I want to include the relative momentum from both particles, but this is a bit hard considering I just assumed 0 momentum from the other party else if(n_macro_cell_particle.get() != (i.get() + 1u) and n_macro_cell_particle.get() > 0u){ // Generally compare } @@ -721,7 +721,7 @@ saw::error_or real_main(int argc, char** argv){ auto& macro_cell = macros.at(index); auto& rho = macro_cell.template get<"pressure">(); - auto& vel = macro_cell.template get<"velocity">(); + auto& vel = macro_cell.template get<"momentum">(); auto& force = macro_cell.template get<"force">(); auto& part_mask = macro_cell.template get<"particle">(); diff --git a/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix b/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix new file mode 100644 index 0000000..d5361c2 --- /dev/null +++ b/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix @@ -0,0 +1,44 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, python3 +, pname +, version +, adaptive-cpp +, kel +, particle_coupling ? "fplbm" +}: + +let + slip_txt = if slip then "slip" else "noslip"; +in stdenv.mkDerivation { + pname = "${pname}-examples-stdkr_2d_gpu_${particle_coupling}"; + 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/schaefer_turek_durst_krause_rannbacher/SConscript b/examples/schaefer_turek_durst_krause_rannbacher/SConscript new file mode 100644 index 0000000..4483d58 --- /dev/null +++ b/examples/schaefer_turek_durst_krause_rannbacher/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/poiseulle_particles_2d_hlbm_gpu', [examples_objects]); + +# Set Alias +env.examples = [ + examples_env.poiseulle_2d_gpu +]; +env.Alias('examples', env.examples); +env.targets += ['examples']; +env.Install('$prefix/bin/', env.examples); diff --git a/examples/schaefer_turek_durst_krause_rannbacher/SConstruct b/examples/schaefer_turek_durst_krause_rannbacher/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/schaefer_turek_durst_krause_rannbacher/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/schaefer_turek_durst_krause_rannbacher/common.hpp b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp new file mode 100644 index 0000000..2d51614 --- /dev/null +++ b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp @@ -0,0 +1,67 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace kel { +namespace lbm { + +/** + * length 2.2m + * viscosity 10^-3 + */ + +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, "momentum">, + Member, "density">, + Member, "porosity"> +>; + +template +using ParticleSpheroidGroup = ParticleGroup>; +} + +} +} diff --git a/examples/schaefer_turek_durst_krause_rannbacher/init.hpp b/examples/schaefer_turek_durst_krause_rannbacher/init.hpp new file mode 100644 index 0000000..b7fdeb0 --- /dev/null +++ b/examples/schaefer_turek_durst_krause_rannbacher/init.hpp @@ -0,0 +1,127 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +template +saw::error_or init( + const converter& conv, + saw::data>& fields, + saw::data>& macros, + saw::data>& particles +){ + (void) conv; + + 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<"momentum">(); + auto& por_f = macros.template get<"porosity">(); + + iterator::apply( + [&](auto& index){ + auto& df = df_f.at(index); + auto& rho = rho_f.at(index); + por_f.at(index).at({}) = {1}; + rho.at({}) = {1}; + auto& vel = vel_f.at(index); + auto eq = equilibrium(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims() + ); + + iterator::apply( + [&](auto& index){ + auto& df = df_f.at(index); + auto& rho = rho_f.at(index); + rho.at({}) = {1}; + auto& vel = vel_f.at(index); + if(info_f.at(index).get() == 2u){ + vel.at({{0u}}) = 0.0; + } + auto eq = equilibrium(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims(), + {{1u,1u}} + ); + + iterator::apply( + [&](auto& index){ + saw::data> middle, ind_vec; + middle.at({{0u}}) = dim_x * 0.25; + middle.at({{1u}}) = dim_y * 0.5; + + ind_vec.at({{0u}}) = index.at({{0u}}).template cast_to(); + ind_vec.at({{1u}}) = index.at({{1u}}).template cast_to(); + + auto dist = middle - ind_vec; + auto dist_2 = saw::math::dot(dist,dist); + if(dist_2.at({}).get() < dim_y*dim_y*0.01){ + porous_f.at(index).at({}) = 0.0; + } + }, + {},// 0-index + df_f.get_dims() + ); + + { + saw::data> radius_p; + radius_p.at({}).set(2); + saw::data> dense_p; + dense_p.at({}).set(1); + particles = create_spheroid_particle_group(radius_p, dense_p, {{16u}}); + } + + return saw::make_void(); +} + +} +} diff --git a/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp new file mode 100644 index 0000000..c7f1014 --- /dev/null +++ b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp @@ -0,0 +1,226 @@ +#include "common.hpp" +#include "init.hpp" +#include "step.hpp" + +#include + +/** + * For deciding what parameters to use maybe? + */ +namespace args { +using namespace saw::schema; + +using LbmArgsStruct = Struct< + Member +>; + +using LbmArgs = Args< + LbmArgsStruct, + Tuple<> +>; +} + +template +saw::error_or lbm_main(const saw::data& args){ + using namespace kel::lbm; + + auto& an = args.template get<"args">(); + + 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 / "stdkr_2d_gpu" / an.template get<"coupling">().stl_view(); + + { + 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_particle_data_ptr = saw::heap>>(); + + std::cout<<"Estimated Bytes: "<,sch::MacroStruct>().get()<(conv,*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_vtk_file(out_dir,"initial_state",0u,*lbm_data_ptr); + if(eov.is_error()){ + return eov; + } + } + + saw::data, encode::Sycl> lbm_sycl_data{sycl_q}; + saw::data, encode::Sycl> lbm_sycl_macro_data{sycl_q}; + saw::data, encode::Sycl> lbm_sycl_particle_data{sycl_q}; + 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); + auto lsdp_view = make_view(lbm_sycl_particle_data); + + 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); + if(eov.is_error()){ + return eov; + } + } + + sycl_q.wait(); + if(i.get() % 64u == 0u){ + { + auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_csv_file(out_dir,"m",i.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + } + // Stream + sycl_q.submit([&](acpp::sycl::handler& h){ + component> stream; + + h.parallel_for(acpp::sycl::range{dim_x,dim_y}, [=](acpp::sycl::id idx){ + saw::data> index; + for(uint64_t i = 0u; i < Desc::D; ++i){ + index.at({{i}}).set(idx[i]); + } + + auto info = info_f.at(index).get(); + + if(info > 0u){ + stream.apply(lsd_view,index,i); + } + }); + }).wait(); + wait.poll(); + + if(print_status){ + std::cout<<"Status: "<().get() * 100 / time_steps.get())<<"%"< 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.stl_view() == "hlbm"){ + return lbm_main(args); + } + + return saw::make_error("Invalid coupling"); +} + + +int main(int argc, char** argv){ + auto eov = k_main(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + std::cerr<<"[Error] "< 0u){ + std::cerr<<" - "< +saw::error_or step( + const converter& conv, + 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 = 16u; + 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(); + + return saw::make_void(); +} + +} +} diff --git a/examples/settling_cubes_2d_ibm_gpu/sim.cpp b/examples/settling_cubes_2d_ibm_gpu/sim.cpp index 33712b5..d7c51bc 100644 --- a/examples/settling_cubes_2d_ibm_gpu/sim.cpp +++ b/examples/settling_cubes_2d_ibm_gpu/sim.cpp @@ -45,7 +45,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "force"> >; @@ -93,7 +93,7 @@ saw::error_or setup_initial_conditions( // auto& df_f = fields.template gStarted hearing about similar cases not long after, many of which were successful on the part of the crooks. Scary stuff. et<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); iterator::apply( [&](auto& index){ @@ -163,7 +163,7 @@ saw::error_or step( auto& parts = p_group.template get<"particles">(); auto& p_mask = p_group.template get<"mask">(); - auto& vels = macros.template get<"velocity">(); + auto& vels = macros.template get<"momentum">(); auto& forces = macros.template get<"force">(); auto p_meta = parts.meta(); @@ -214,7 +214,7 @@ saw::error_or step( // Lagrange indicator position auto p_pos_lag = p_pos + transformed_pos; - // Pick the closest velocity and clamp it + // Pick the closest momentum and clamp it saw::data> p_cell_pos; saw::data> p_cell_pos_vec; for(uint64_t i{0u}; i < Desc::D; ++i){ diff --git a/examples/settling_spheres_2d_hlbm_gpu/sim.cpp b/examples/settling_spheres_2d_hlbm_gpu/sim.cpp index f7e47b6..5913a1f 100644 --- a/examples/settling_spheres_2d_hlbm_gpu/sim.cpp +++ b/examples/settling_spheres_2d_hlbm_gpu/sim.cpp @@ -49,7 +49,7 @@ using PorChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "porosity"> >; @@ -91,7 +91,7 @@ saw::error_or setup_initial_conditions( // auto& df_f = fields.template gStarted hearing about similar cases not long after, many of which were successful on the part of the crooks. Scary stuff. et<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); iterator::apply( [&](auto& index){ @@ -168,7 +168,7 @@ saw::error_or step( auto& parts = p_group.template get<"particles">(); auto& p_mask = p_group.template get<"mask">(); - auto& vels = macros.template get<"velocity">(); + auto& vels = macros.template get<"momentum">(); auto& forces = macros.template get<"force">(); auto p_meta = parts.meta(); diff --git a/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp b/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp index 9d2ec12..bbe8cfd 100644 --- a/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp +++ b/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp @@ -49,7 +49,7 @@ using ForceChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density"> >; @@ -122,7 +122,7 @@ saw::error_or setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); iterator::apply( [&](auto& index){ diff --git a/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp index 9e96219..9cf974a 100644 --- a/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp +++ b/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp @@ -49,7 +49,7 @@ using ForceChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "force">, Member, "porosity"> @@ -125,7 +125,7 @@ saw::error_or setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator::apply( diff --git a/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp index 9c28d88..b63a003 100644 --- a/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp +++ b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp @@ -46,7 +46,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "force">, Member, "porosity"> @@ -122,7 +122,7 @@ saw::error_or setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator::apply( diff --git a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp index 75cbe54..6b86c68 100644 --- a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp +++ b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp @@ -47,7 +47,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density">, Member, "force">, Member, "porosity"> @@ -124,7 +124,7 @@ saw::error_or setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator::apply( diff --git a/modules/core/c++/boundary.hpp b/modules/core/c++/boundary.hpp index 7f185ad..e7567c9 100644 --- a/modules/core/c++/boundary.hpp +++ b/modules/core/c++/boundary.hpp @@ -133,14 +133,14 @@ template class component final { private: saw::data> density_; - saw::data> velocity_; + saw::data> momentum_; public: component( saw::data> density__, - saw::data> velocity__ + saw::data> momentum__ ): density_{density__}, - velocity_{velocity__} + momentum_{momentum__} {} template @@ -151,7 +151,7 @@ public: auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); - auto eq = equilibrium(density_,velocity_); + auto eq = equilibrium(density_,momentum_); dfs_old_f.at(index) = eq; } @@ -726,12 +726,12 @@ public: template class component, Encode> final { private: - saw::data> velocity_; + saw::data> momentum_; public: component( - saw::data> velocity__ + saw::data> momentum__ ): - velocity_{velocity__} + momentum_{momentum__} {} template @@ -753,20 +753,20 @@ public: } saw::data> rho; - rho.at({}) = (dfs.at({0u}) + dfs.at({4u}) + dfs.at({3u}) + dir_sum * 2) / (velocity_.at({{0u}}) + saw::data{static_cast::type>(1.0)}); + rho.at({}) = (dfs.at({0u}) + dfs.at({4u}) + dfs.at({3u}) + dir_sum * 2) / (momentum_.at({{0u}}) + saw::data{static_cast::type>(1.0)}); if constexpr (East) { - dfs.at({2u}) = dfs.at({1u}) + saw::data{static_cast::type>(2.0 / 3.0)} * rho.at({}) * velocity_.at({{0u}}); - dfs.at({6u}) = dfs.at({5u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}}) - + (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * velocity_.at({{1u}})) * 0.5f; - dfs.at({8u}) = dfs.at({7u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}}) - + (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * velocity_.at({{1u}})) * 0.5f; + dfs.at({2u}) = dfs.at({1u}) + saw::data{static_cast::type>(2.0 / 3.0)} * rho.at({}) * momentum_.at({{0u}}); + dfs.at({6u}) = dfs.at({5u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}}) + + (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * momentum_.at({{1u}})) * 0.5f; + dfs.at({8u}) = dfs.at({7u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}}) + + (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * momentum_.at({{1u}})) * 0.5f; }else if constexpr ( not East ){ - dfs.at({1u}) = dfs.at({2u}) + saw::data{static_cast::type>(2.0 / 3.0)} * rho.at({}) * velocity_.at({{0u}}); - dfs.at({5u}) = dfs.at({6u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}}) - + (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * velocity_.at({{1u}}))*0.5f; - dfs.at({7u}) = dfs.at({8u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}}) - + (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * velocity_.at({{1u}}))*0.5f; + dfs.at({1u}) = dfs.at({2u}) + saw::data{static_cast::type>(2.0 / 3.0)} * rho.at({}) * momentum_.at({{0u}}); + dfs.at({5u}) = dfs.at({6u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}}) + + (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * momentum_.at({{1u}}))*0.5f; + dfs.at({7u}) = dfs.at({8u}) + saw::data{static_cast::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}}) + + (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * momentum_.at({{1u}}))*0.5f; } } }; diff --git a/modules/core/c++/collision.hpp b/modules/core/c++/collision.hpp index 023f61f..39514b8 100644 --- a/modules/core/c++/collision.hpp +++ b/modules/core/c++/collision.hpp @@ -75,7 +75,7 @@ public: auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data>& rho = rho_f.at(index); saw::data>& vel = vel_f.at(index); @@ -135,7 +135,7 @@ public: auto& dfs = dfs_old_f.at(index); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data>& rho = rho_f.at(index); diff --git a/modules/core/c++/component/psm.hpp b/modules/core/c++/component/psm.hpp index cada53c..5832760 100644 --- a/modules/core/c++/component/psm.hpp +++ b/modules/core/c++/component/psm.hpp @@ -39,7 +39,7 @@ public: auto& porous_f = macros.template get<"porosity">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data>& rho = rho_f.at(index); saw::data>& vel = vel_f.at(index); diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp index e2527ce..382d067 100644 --- a/modules/core/c++/fplbm.hpp +++ b/modules/core/c++/fplbm.hpp @@ -46,7 +46,7 @@ public: auto& rho_f = macros.template get<"density">(); auto& rho = rho_f.at(index); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& vel = vel_f.at(index); compute_rho_u(dfs,rho,vel); @@ -86,7 +86,7 @@ public: auto& rho_f = macros.template get<"density">(); saw::data>& rho = rho_f.at(index); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data> vel = vel_f.at(index); auto& force_f = macros.template get<"force">(); @@ -148,7 +148,7 @@ public: // auto& dfs = dfs_old_f.at(index); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); auto& force_f = macros.template get<"force">(); @@ -245,7 +245,7 @@ public: auto& dfs = dfs_old_f.at(index); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); // Temporary find a better way diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp index cc21dcf..e8ea66a 100644 --- a/modules/core/c++/hlbm.hpp +++ b/modules/core/c++/hlbm.hpp @@ -70,7 +70,7 @@ public: auto& porosity_f = macros.template get<"porosity">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data>& rho = rho_f.at(index); saw::data>& vel = vel_f.at(index); @@ -125,7 +125,7 @@ public: auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); auto& part_spheroid_group = part_group; - auto& mvel_f = macros.template get<"velocity">(); + auto& mvel_f = macros.template get<"momentum">(); auto& mrho_f = macros.template get<"density">(); auto& mpor_f = macros.template get<"porosity">(); diff --git a/modules/core/c++/macroscopic.hpp b/modules/core/c++/macroscopic.hpp index 19ab3e9..768fef0 100644 --- a/modules/core/c++/macroscopic.hpp +++ b/modules/core/c++/macroscopic.hpp @@ -28,6 +28,7 @@ void compute_rho_u ( } } + // TODO - Change to momentum, but that needs checks in every impacted function for(size_t i = 0; i < Desc::D; ++i){ vel().at({{i}}) = vel().at({{i}}) / rho().at({}); } diff --git a/modules/core/c++/particle.hpp b/modules/core/c++/particle.hpp index 2c60884..efcc3fb 100644 --- a/modules/core/c++/particle.hpp +++ b/modules/core/c++/particle.hpp @@ -25,7 +25,7 @@ public: // Compute forces - // Update particle velocity + // Update particle momentum verlet_step_lambda_old(p,{1.0}); // Update porosity over lattice nodes diff --git a/modules/core/c++/psm.hpp b/modules/core/c++/psm.hpp index 15ad0f8..2daff5e 100644 --- a/modules/core/c++/psm.hpp +++ b/modules/core/c++/psm.hpp @@ -55,7 +55,7 @@ public: auto& porous_f = macros.template get<"porosity">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data>& rho = rho_f.at(index); saw::data>& vel = vel_f.at(index); @@ -105,116 +105,418 @@ public: component() = default; template - void apply(const saw::data& field, const saw::data& macros, const saw::data& pg, saw::data> index, saw::data time_step, saw::data sub_steps) const { - - using dfi = df_info; - bool is_even = ((time_step.get() % 2) == 0); - - saw::data> one; - one.at({}) = 1.0; +void apply( + const saw::data& field, + const saw::data& macros, + const saw::data& pg, + saw::data> index, + saw::data time_step, + saw::data sub_steps +) const { - auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + using dfi = df_info; - auto& por_f = macros.template get<"porosity">(); - auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); - auto& force_f = macros.template get<"force">(); - - { - auto parts = pg.template get<"particles">(); - auto parts_size = parts.meta().at({0u}); - - auto& p_coll = pg.template get<"collision">().at({}); - auto& p_rad = p_coll.template get<"radius">(); - - auto& pi = parts.at(index); - auto& pirb = pi.template get<"rigid_body">(); - saw::data>& pirb_pos = pirb.template get<"position">(); - auto& pirb_pos_old = pirb.template get<"position_old">(); - - saw::data> ts; - ts.at({}) = one.at({}) / sub_steps.template cast_to(); - - saw::data> start; - saw::data> stop; - - auto eo_aabb = particle_aabb::calculate(pg,{{0u}},vel_f.meta()); - if(eo_aabb.is_error()){ - return; - } - auto& aabb = eo_aabb.get_value(); + bool is_even = ((time_step.get() % 2u) == 0u); - /// Ok, I iterate over the space which covers our particle? So lower bounds to upper bounds - start = aabb.template get<"a">(); - stop = aabb.template get<"b">(); + /* + * ------------------------------------------------------------ + * Fields + * ------------------------------------------------------------ + */ - saw::data> force_p; - for(uint64_t i{0u}; i < Descriptor::D; ++i){ - force_p.at({{i}}) = 0.0; - } - auto vel_p_old = (pirb_pos-pirb_pos_old) / ts; - - iterator::apply([&](const auto& index_f) -> void{ - // ask for the d_k value here. - // For every value im iterating over I need sth - // std::cout<<"Pos: "<> rel_dist = saw::math::vectorize_data(index_f).template cast_to() - pirb_pos; - saw::data> eps; - eps.at({}) = 1.5f; - - auto& por = por_f.at(index_f); - por = particle_porosity>::calculate(rel_dist,p_rad,eps); - - if(por.at({}).get() >= 1.0f){ - return; - } - - saw::data> momentum; - for(uint64_t i{0u}; i < Descriptor::D; ++i){ - momentum.at({{i}}) = 0.0; - } - - for(uint64_t i{0u}; i < Descriptor::Q; ++i){ - saw::data> e_i; - saw::data> n_ind_i; - for(uint64_t k{0u}; k < Descriptor::D; ++k){ - e_i.at({{k}}) = (dfi::directions[i])[k]; - n_ind_i.at({k}) = index_f.at({k}) + (dfi::directions[i])[k]; - } - - uint64_t i_opp = dfi::opposite_index[i]; - - auto u_p_e = saw::math::dot(e_i,vel_p_old); - - saw::data dfs_added = dfs.at({i})*(saw::data{1}-u_p_e.at({})) + dfs_old_f.at(n_ind_i).at({i_opp})*(saw::data{1}+u_p_e.at({})); - saw::data> dfs_added_v; - dfs_added_v.at({}) = dfs_added; - auto ei_dfs = e_i * dfs_added_v; - - momentum = momentum + ei_dfs; - } - // technically needs to adjust for rotation as well - - auto& force = force_f.at(index_f); - auto& rho = rho_f.at(index_f); - // To Fluid - - auto flip_por = one - por; - force = momentum * flip_por; - // To Particle - force_p = force_p - force; - },start,stop); - - auto& pirb_acc = pirb.template get<"acceleration">(); - pirb_acc = force_p / pg.template get<"total_mass">().at({}); - - for(saw::data i{0u}; i < sub_steps; ++i){ - verlet_step_lambda(pi,ts); - } - } - } + auto& dfs_old_f = + (is_even) + ? field.template get<"dfs_old">() + : field.template get<"dfs">(); + + auto& por_f = macros.template get<"porosity">(); + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"momentum">(); + auto& force_f = macros.template get<"force">(); + + saw::data> one; + one.at({}) = 1.0; + + saw::data> half; + half.at({}) = 0.5; + + /* + * ------------------------------------------------------------ + * Particle + * ------------------------------------------------------------ + */ + + { + auto parts = pg.template get<"particles">(); + + auto& p_coll = + pg.template get<"collision">().at({}); + + auto& p_rad = + p_coll.template get<"radius">(); + + auto& pi = parts.at(index); + + auto& pirb = + pi.template get<"rigid_body">(); + + saw::data>& pirb_pos = + pirb.template get<"position">(); + + auto& pirb_pos_old = + pirb.template get<"position_old">(); + + + /* + * -------------------------------------------------------- + * Particle timestep + * -------------------------------------------------------- + */ + + saw::data> ts; + + ts.at({}) = + one.at({}) + / sub_steps.template cast_to(); + + + /* + * -------------------------------------------------------- + * Particle momentum + * -------------------------------------------------------- + */ + + auto vel_p_old = + (pirb_pos - pirb_pos_old) / ts; + + + /* + * -------------------------------------------------------- + * Find particle AABB + * -------------------------------------------------------- + */ + + auto eo_aabb = + particle_aabb< + typename ParticleSchema::ValueType + >::calculate( + pg, + {{0u}}, + vel_f.meta() + ); + + if(eo_aabb.is_error()){ + return; + } + + auto& aabb = eo_aabb.get_value(); + + saw::data> start; + saw::data> stop; + + start = aabb.template get<"a">(); + stop = aabb.template get<"b">(); + + + /* + * -------------------------------------------------------- + * Total force acting on particle + * -------------------------------------------------------- + */ + + saw::data> force_p; + + for(uint64_t k = 0u; k < Descriptor::D; ++k){ + force_p.at({{k}}) = 0.0; + } + + + /* + * -------------------------------------------------------- + * Iterate over particle AABB + * -------------------------------------------------------- + */ + + iterator::apply( + [&](const auto& index_f) -> void { + + /* + * ------------------------------------------------ + * Distribution function at current cell + * ------------------------------------------------ + */ + + auto& dfs = + dfs_old_f.at(index_f); + + + /* + * ------------------------------------------------ + * Particle relative position + * ------------------------------------------------ + */ + + saw::data> rel_dist = + saw::math::vectorize_data(index_f) + .template cast_to() + - pirb_pos; + + + /* + * ------------------------------------------------ + * PSM particle porosity + * ------------------------------------------------ + */ + + saw::data> eps; + eps.at({}) = 1.5; + + auto& por = + por_f.at(index_f); + + por = + particle_porosity< + T, + Descriptor::D, + 1u, + por::ParticleSpheroid + >::calculate( + rel_dist, + p_rad, + eps + ); + + + /* + * Outside particle. + */ + + if(por.at({}).get() >= 1.0){ + return; + } + + + /* + * ------------------------------------------------ + * Momentum exchange for this cell + * ------------------------------------------------ + */ + + saw::data> momentum; + + for(uint64_t k = 0u; k < Descriptor::D; ++k){ + momentum.at({{k}}) = 0.0; + } + + + /* + * ------------------------------------------------ + * D3Q27 / general Q momentum exchange + * + * We exclude i=0 because the rest population has + * no momentum. + * + * The original code summed both directions of + * every opposite pair. Therefore we apply 1/2 + * after the sum. + * ------------------------------------------------ + */ + + for(uint64_t i = 1u; i < Descriptor::Q; ++i){ + + uint64_t i_opp = + dfi::opposite_index[i]; + + + /* + * Direction vector e_i. + */ + + saw::data> e_i; + + saw::data< + sch::FixedArray< + sch::UInt64, + Descriptor::D + > + > n_ind_i; + + + for(uint64_t k = 0u; + k < Descriptor::D; + ++k){ + + e_i.at({{k}}) = + dfi::directions[i][k]; + + /* + * Neighbor in direction e_i. + */ + + n_ind_i.at({k}) = + index_f.at({k}) + + dfi::directions[i][k]; + } + + + /* + * ------------------------------------------------ + * Particle momentum projected onto e_i + * ------------------------------------------------ + */ + + auto u_p_e = + saw::math::dot( + e_i, + vel_p_old + ); + + + /* + * ------------------------------------------------ + * Momentum-exchange population + * + * f_i(x) + * + * + f_-i(x + e_i) + * + * with moving-wall correction. + * ------------------------------------------------ + */ + + saw::data dfs_added = + dfs.at({i}) + * ( + saw::data{1.0} + - u_p_e.at({}) + ) + + + dfs_old_f + .at(n_ind_i) + .at({i_opp}) + * ( + saw::data{1.0} + + u_p_e.at({}) + ); + + + /* + * Convert scalar into vector multiplication. + */ + + saw::data> dfs_added_v; + + dfs_added_v.at({}) = + dfs_added; + + + /* + * ------------------------------------------------ + * Delta p_i + * + * ------------------------------------------------ + */ + + momentum = + momentum + + e_i * dfs_added_v; + } + + + /* + * ------------------------------------------------ + * Since both i and i_opp were included above, + * compensate for double counting. + * + * This is the important change from your original + * implementation. + * ------------------------------------------------ + */ + + + + momentum = + momentum * half; + + + /* + * ------------------------------------------------ + * PSM solid fraction + * + * por = 1 -> fluid + * por = 0 -> particle + * + * Therefore: + * + * flip_por = 1 - por + * ------------------------------------------------ + */ + + auto flip_por = + one - por; + + + /* + * ------------------------------------------------ + * Force transferred TO FLUID + * ------------------------------------------------ + */ + + auto& force = + force_f.at(index_f); + + force = + momentum * flip_por; + + + /* + * ------------------------------------------------ + * Force transferred TO PARTICLE + * + * Newton's third law. + * ------------------------------------------------ + */ + + force_p = + force_p - force; + + }, + start, + stop + ); + + + /* + * ------------------------------------------------------------ + * Particle acceleration + * ------------------------------------------------------------ + */ + + auto& pirb_acc = + pirb.template get<"acceleration">(); + + pirb_acc = + force_p + / pg.template get<"total_mass">().at({}); + + + /* + * ------------------------------------------------------------ + * Integrate particle + * ------------------------------------------------------------ + */ + + for(saw::data i{0u}; + i < sub_steps; + ++i){ + + verlet_step_lambda( + pi, + ts + ); + } + } +} }; } diff --git a/modules/core/c++/write_csv.hpp b/modules/core/c++/write_csv.hpp index a60e208..aa28613 100644 --- a/modules/core/c++/write_csv.hpp +++ b/modules/core/c++/write_csv.hpp @@ -91,7 +91,7 @@ struct lbm_csv_writer> { template struct lbm_csv_writer> { static saw::error_or apply(std::ostream& csv_file, const saw::data>& field){ - static_assert(D > 0, "Non-dimensionality is bad for velocity."); + static_assert(D > 0, "Non-dimensionality is bad for momentum."); // csv_file<<"VECTORS "<> { return saw::make_void(); } static saw::error_or apply(std::ostream& vtk_file, const saw::data>& field){ - static_assert(D > 0, "Non-dimensionality is bad for velocity."); + static_assert(D > 0, "Non-dimensionality is bad for momentum."); static_assert(D <= 3, "4th dimension as well. Mostly due to vtk."); // vtk_file<<"VECTORS "<(rho,vel); for(saw::data i{0u}; i.get() < Descriptor::Q; ++i){ - SAW_EXPECT(eq.at(i).get() == dfi::weights[i.get()], std::string{"No velocity and normalized rho should be exactly the weights: "} + std::to_string(eq.at(i).get()) + std::string{" "} + std::to_string(dfi::weights[i.get()])); + SAW_EXPECT(eq.at(i).get() == dfi::weights[i.get()], std::string{"No momentum and normalized rho should be exactly the weights: "} + std::to_string(eq.at(i).get()) + std::string{" "} + std::to_string(dfi::weights[i.get()])); } } diff --git a/modules/core/tests/vtk_write.cpp b/modules/core/tests/vtk_write.cpp index c99d752..1de25bf 100644 --- a/modules/core/tests/vtk_write.cpp +++ b/modules/core/tests/vtk_write.cpp @@ -31,7 +31,7 @@ using CellStruct = Struct< template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member >; @@ -47,7 +47,7 @@ SAW_TEST("VTK Write test example"){ saw::data, 2>> cells{{{2u,2u}}}; auto& cell_0 = cells.at({{{0,0}}}); - cell_0.template get<"velocity">()= {{0.5,-0.1}}; + cell_0.template get<"momentum">()= {{0.5,-0.1}}; cell_0.template get<"pressure">().set(1.1); auto eov = lbm::impl::lbm_vtk_writer, 2>>::apply(sstream, cells); @@ -56,7 +56,7 @@ SAW_TEST("VTK Write test example"){ // I want to print it to see it for myself. For now I have no tooling to more easily view associated and potentially generated files std::cout<::type; @@ -74,7 +74,7 @@ SAW_TEST("VTK Write raw test example"){ auto meta = cells.dims(); auto& cell_0 = cells.at({{{0,0}}}); - cell_0.template get<"velocity">()= {{0.5,-0.1}}; + cell_0.template get<"momentum">()= {{0.5,-0.1}}; cell_0.template get<"pressure">().set(1.1); auto eov = lbm::impl::lbm_vtk_writer_raw,2u>::apply(sstream, &cell_0, meta); @@ -83,7 +83,7 @@ SAW_TEST("VTK Write raw test example"){ // I want to print it to see it for myself. For now I have no tooling to more easily view associated and potentially generated files std::cout<::type; diff --git a/poc/triangle_stream/common.hpp b/poc/triangle_stream/common.hpp index 927f8f8..848d337 100644 --- a/poc/triangle_stream/common.hpp +++ b/poc/triangle_stream/common.hpp @@ -37,7 +37,7 @@ using RhoChunk = Chunk, 0u, dim_x, dim_y>; template using MacroStruct = Struct< - Member, "velocity">, + Member, "momentum">, Member, "density"> >; } -- cgit v1.2.3 From c634aac438972b7bafb5fbbeba7962d174cc0d03 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Fri, 14 Aug 2026 23:43:36 +0200 Subject: Experimental validation --- default.nix | 5 + .../.nix/derivation.nix | 1 - .../SConscript | 2 +- .../common.hpp | 7 +- .../init.hpp | 2 +- .../schaefer_turek_durst_krause_rannbacher/sim.cpp | 16 ++- .../step.hpp | 125 +++++++++++---------- 7 files changed, 86 insertions(+), 72 deletions(-) diff --git a/default.nix b/default.nix index f560d7e..d7aeb9c 100644 --- a/default.nix +++ b/default.nix @@ -218,6 +218,11 @@ in rec { inherit kel; }; + stdkr_2d = pkgs.callPackage ./examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix { + inherit pname version stdenv forstio adaptive-cpp; + inherit kel; + }; + heterogeneous_computing = pkgs.callPackage ./examples/heterogeneous_computing/.nix/derivation.nix { inherit pname version stdenv forstio adaptive-cpp; inherit kel; diff --git a/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix b/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix index d5361c2..958955f 100644 --- a/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix +++ b/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix @@ -12,7 +12,6 @@ }: let - slip_txt = if slip then "slip" else "noslip"; in stdenv.mkDerivation { pname = "${pname}-examples-stdkr_2d_gpu_${particle_coupling}"; inherit version; diff --git a/examples/schaefer_turek_durst_krause_rannbacher/SConscript b/examples/schaefer_turek_durst_krause_rannbacher/SConscript index 4483d58..b33b7ae 100644 --- a/examples/schaefer_turek_durst_krause_rannbacher/SConscript +++ b/examples/schaefer_turek_durst_krause_rannbacher/SConscript @@ -23,7 +23,7 @@ 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/poiseulle_particles_2d_hlbm_gpu', [examples_objects]); +examples_env.poiseulle_2d_gpu = examples_env.Program('#bin/stdkr_2d', [examples_objects]); # Set Alias env.examples = [ diff --git a/examples/schaefer_turek_durst_krause_rannbacher/common.hpp b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp index 2d51614..ec6f80f 100644 --- a/examples/schaefer_turek_durst_krause_rannbacher/common.hpp +++ b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp @@ -15,11 +15,14 @@ namespace lbm { /** * length 2.2m + * height 0.41m * viscosity 10^-3 + * particle diameter 0.1m + * position at (0.16m,0.2m) */ -constexpr uint64_t dim_y = 256ul; -constexpr uint64_t dim_x = dim_y * 20ul; +constexpr uint64_t dim_y = 512ul; +constexpr uint64_t dim_x = static_cast(dim_y * (220.0 / 41.0)); constexpr uint64_t particle_amount = 1ul; diff --git a/examples/schaefer_turek_durst_krause_rannbacher/init.hpp b/examples/schaefer_turek_durst_krause_rannbacher/init.hpp index b7fdeb0..b0e4020 100644 --- a/examples/schaefer_turek_durst_krause_rannbacher/init.hpp +++ b/examples/schaefer_turek_durst_krause_rannbacher/init.hpp @@ -114,7 +114,7 @@ saw::error_or init( { saw::data> radius_p; - radius_p.at({}).set(2); + radius_p.at({}).set(conv.meter_si_to_lbm({{0.16}}).handle().get()); saw::data> dense_p; dense_p.at({}).set(1); particles = create_spheroid_particle_group(radius_p, dense_p, {{16u}}); diff --git a/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp index c7f1014..6113f6c 100644 --- a/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp +++ b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp @@ -46,12 +46,12 @@ saw::error_or lbm_main(const saw::data& args){ converter conv { // delta_x - {{1.0}}, + {{0.41 / dim_y}}, // delta_t - {{1.0}} + {{1.067e-4}} }; - print_lbm_meta(conv,{0.1},{1e-4},{0.4 * dim_y}); + print_lbm_meta(conv,{1e-3},{1e-4},{0.5 * dim_y}); // saw::data> meta{{dim_x,dim_y}}; auto lbm_data_ptr = saw::heap>>(); @@ -138,7 +138,7 @@ saw::error_or lbm_main(const saw::data& args){ } } { - auto eov = write_csv_file(out_dir,"m",i.get(), *lbm_macro_data_ptr); + auto eov = write_csv_file(out_dir,"stdkr_2d",i.get(), *lbm_macro_data_ptr); if(eov.is_error()){ return eov; } @@ -178,7 +178,7 @@ saw::error_or lbm_main(const saw::data& args){ } } { - auto eov = write_vtk_file(out_dir,"m",time_steps.get(), *lbm_macro_data_ptr); + auto eov = write_vtk_file(out_dir,"stdkr_2d",time_steps.get(), *lbm_macro_data_ptr); if(eov.is_error()){ return eov; } @@ -197,13 +197,17 @@ saw::error_or k_main(int argc, char** argv){ saw::data args; { auto eov = saw::parse_args(args,argc,argv); - return eov; + if(eov.is_error()){ + return eov; + } } auto& an = args.template get<"args">(); auto& coupling = an.template get<"coupling">(); if(coupling.stl_view() == "hlbm"){ return lbm_main(args); + } else if(coupling.stl_view() == "fplbm"){ + return lbm_main(args); } return saw::make_error("Invalid coupling"); diff --git a/examples/schaefer_turek_durst_krause_rannbacher/step.hpp b/examples/schaefer_turek_durst_krause_rannbacher/step.hpp index 37eb4e7..1a70f6a 100644 --- a/examples/schaefer_turek_durst_krause_rannbacher/step.hpp +++ b/examples/schaefer_turek_durst_krause_rannbacher/step.hpp @@ -14,81 +14,84 @@ saw::error_or step( saw::data t_i, device& dev ){ + static constexpr saw::data sub_steps{0u}; + 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]); - } + if constexpr ( std::is_same_v ){ + 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(); + hlbm_reset.apply(fields,macros,index,t_i); + }); + }).wait(); + q.submit([&](acpp::sycl::handler& h){ + component> hlbm_one_part; - 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]); + } - 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,{0u}); + }); + }).wait(); - hlbm_one_part.apply(fields,macros,particles,index,t_i,{16u}); - }); - }).wait(); + q.submit([&](acpp::sycl::handler& h){ + component> collision{1.0}; + component> bb; - // 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 = 16u; + if(t_i.get() < target_t_i){ + return 1.0 + (0.01 / target_t_i) * t_i.get(); + } + return 1.001; + }() + }; + component,encode::Sycl> flow_out{1.0}; - component,encode::Sycl> flow_in{ - [&](){ - uint64_t target_t_i = 16u; - if(t_i.get() < target_t_i){ - return 1.0 + (0.01 / target_t_i) * t_i.get(); + 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]); } - 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); + auto info = info_f.at(index); + + switch(info.get()){ + case 0u: break; - case 3u: - flow_in.apply(fields,index,t_i); - collision.apply(fields,macros,index,t_i); + case 1u: + bb.apply(fields,index,t_i); break; - case 4u: - flow_out.apply(fields,index,t_i); - collision.apply(fields,macros,index,t_i); - break; - default: - break; - } - }); - }).wait(); + 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(); + }else if constexpr ( std::is_same_v ){ + } return saw::make_void(); } -- cgit v1.2.3