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 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 55 insertions(+), 3 deletions(-) (limited to 'modules/core/c++/boundary.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) { -- 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 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 477 insertions(+), 28 deletions(-) (limited to 'modules/core/c++/boundary.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); + } } - } }; -- 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 --- modules/core/c++/boundary.hpp | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) (limited to 'modules/core/c++/boundary.hpp') 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; } } }; -- cgit v1.2.3