From d25750adaeec754fbad0e79a5570ed3b5dc02513 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Thu, 17 Apr 2025 18:06:54 +0200 Subject: Moved examples and added some fun rendering --- SConstruct | 6 + c++/SConscript | 18 +- c++/examples/cavity_2d.cpp | 531 ------------------------------------------ examples/SConscript | 33 +++ examples/cavity_2d.cpp | 564 +++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 605 insertions(+), 547 deletions(-) delete mode 100644 c++/examples/cavity_2d.cpp create mode 100644 examples/SConscript create mode 100644 examples/cavity_2d.cpp diff --git a/SConstruct b/SConstruct index 71818ae..32bf874 100644 --- a/SConstruct +++ b/SConstruct @@ -43,6 +43,11 @@ env_vars.Add('prefix', 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'], @@ -58,6 +63,7 @@ env.targets = []; Export('env') SConscript('c++/SConscript') +SConscript('examples/SConscript') SConscript('tests/SConscript') env.Alias('cdb', env.cdb); diff --git a/c++/SConscript b/c++/SConscript index 772d526..57091af 100644 --- a/c++/SConscript +++ b/c++/SConscript @@ -19,21 +19,7 @@ env.sources += core_env.sources; env.headers += core_env.headers; -## Shared lib +## Static lib objects = [] core_env.add_source_files(objects, core_env.sources, shared=False); - -core_env.examples = []; -# Cavity2D -core_env.cavity_2d_source = sorted(glob.glob(dir_path + "/examples/cavity_2d.cpp")); -env.sources += core_env.cavity_2d_source; -core_env.cavity_2d = core_env.Program('#bin/cavity_2d', [core_env.cavity_2d_source, core_env.objects]); -core_env.examples += core_env.cavity_2d; - -# Set Alias -env.Alias('examples', [core_env.examples]); - -env.targets += ['examples']; - -# Install -env.Install('$prefix/bin', [core_env.examples]); +env.library_static = core_env.StaticLibrary('#build/kel-lbm', [objects]); diff --git a/c++/examples/cavity_2d.cpp b/c++/examples/cavity_2d.cpp deleted file mode 100644 index 00c1741..0000000 --- a/c++/examples/cavity_2d.cpp +++ /dev/null @@ -1,531 +0,0 @@ -#include "../descriptor.hpp" - -/** - */ -#include - -#include -#include -#include - -namespace kel { -namespace lbm { -namespace sch { -using namespace saw::schema; - -/** - * Basic distribution function - * Base type - * D - * Q - * Scalar factor - * D factor - * Q factor - */ -using T = Float32; -using D2Q5 = Descriptor<2,5>; -using D2Q9 = Descriptor<2,9>; - -template -using DfCell = Cell; - -template -using CellInfo = Cell; - -/** - * Basic type for simulation - */ -template -using CellStruct = Struct< - Member, "dfs">, - Member, "dfs_old">, - Member, "info"> ->; - - -using CavityFieldD2Q9 = CellField>; -} - -/** - * Calculate the macroscopic variables rho and u in Lattice Units. - */ -template -void compute_rho_u ( - saw::data>& dfs, - typename saw::native_data_type::type& rho, - std::array::type, 2>& vel - ) -{ - using dfi = df_info; - - rho = 0; - std::fill(vel.begin(), vel.end(), 0); - - for(size_t i = 0; i < Desc::Q; ++i){ - rho += dfs(i).get(); - vel[0] += dfi::directions[i][0] * dfs(i).get(); - vel[1] += dfi::directions[i][1] * dfs(i).get(); - } - - vel[0] /= rho; - vel[1] /= rho; -} - -/** - * Unsure if feasible. I would have a rho normalization on the dfs and then I would use the const rho_computation - */ -template -void compute_const_rho_u ( - saw::data>& dfs, - const typename saw::native_data_type::type rho, - std::array::type, 2>& vel - ) -{ - using dfi = df_info; - - saw::native_data_type::type real_rho = 0; - std::fill(vel.begin(), vel.end(), 0); - - for(size_t i = 0; i < Desc::Q; ++i){ - real_rho += dfs(i).get(); - vel[0] += dfi::directions[i][0] * dfs(i).get(); - vel[1] += dfi::directions[i][1] * dfs(i).get(); - } - for(size_t i = 0; i < Desc::Q; ++i){ - dfs(i).set(dfs(i).get() * (rho/real_rho)); - } - - vel[0] *= real_rho / (rho*rho); - vel[1] *= real_rho / (rho*rho); -} - -/** - * Calculates the equilibrium for each direction - */ -template -std::array::type,Desc::Q> equilibrium( - typename saw::native_data_type::type rho, - const std::array::type, Desc::D>& vel -){ - using dfi = df_info; - - typename std::array::type,Desc::Q> eq; - - for(std::size_t i = 0; i < eq.size(); ++i){ - auto vel_c = (vel[0]*dfi::directions[i][0] + vel[1]*dfi::directions[i][1]); - auto vel_c_cs2 = vel_c * dfi::inv_cs2; - eq[i] = dfi::weights[i] * rho * ( - 1.0 - + vel_c_cs2 - - dfi::inv_cs2 * 0.5 * ( vel[0] * vel[0] + vel[1] * vel[1] ) - + vel_c_cs2 * vel_c_cs2 * 0.5 - ); - } - - return eq; -} - -/** - * A reason for why a component based design is really good can be seen in my LR solver example - * - * Add Expression Templates and you're golden. - */ -template -class component; - -/* -template -class df_cell_view; -*/ -/** - * Minor helper for the AA-Pull Pattern, so I can use only one lattice - * - * Am I sure I want to use AA this way? - * Esoteric Twist technically reduces the needed memory access footprint - */ -/* -template -class df_cell_view, Encode> { -public: - using Schema = sch::Cell; -private: - std::array::type>*, QN> view_; -public: - df_cell_view(const std::array::type>*, QN>& view): - view_{view} - {} -}; -*/ -namespace cmpt { -struct BounceBack{}; -struct MovingWall {}; -struct BGK {}; -struct ConstRhoBGK {}; -} - - -/** - * Full-Way BounceBack. I suspect that the moving wall requires half-way bounce back. - */ -template -class component { -public: - - void apply(saw::data>& dfs){ - using dfi = df_info; - - // Technically use .copy() - auto df_cpy = dfs; - - for(uint64_t i = 1u; i < Desc::Q; ++i){ - dfs({i}) = df_cpy({dfi::opposite_index.at(i)}); - } - } -}; - -/** - * Full-Way moving wall Bounce back, something is not right here. - * Technically it should reflect properly. - */ -template -class component { -public: - std::array::type, Desc::D> lid_vel; - -public: - void apply( - saw::data>& dfs - ){ - using dfi = df_info; - - // Technically use .copy() - /* - auto dfs_cpy = dfs; - - for(uint64_t i = 0u; i < Desc::Q; ++i){ - dfs({dfi::opposite_index.at(i)}) = dfs_cpy({i}) - 2.0 * dfi::weights[i] * 1.0 * ( lid_vel[0] * dfi::directions[i][0] + lid_vel[1] * dfi::directions[i][1]) * dfi::inv_cs2; - } - */ - } -}; - -template -class component { -public: - typename saw::native_data_type::type relaxation_; -public: - void apply(saw::data>& dfs){ - typename saw::native_data_type::type rho; - std::array::type, Desc::D> vel; - compute_rho_u(dfs,rho,vel); - auto eq = equilibrium(rho,vel); - - for(uint64_t i = 0u; i < Desc::Q; ++i){ - dfs({i}).set(dfs({i}).get() + (1.0 / relaxation_) * (eq[i] - dfs({i}).get())); - } - } -}; - -template -class component { -public: - typename saw::native_data_type::type relaxation_; - typename saw::native_data_type::type rho_; -public: - void apply(saw::data>& dfs){ - std::array::type, Desc::D> vel; - compute_const_rho_u(dfs,rho_,vel); - auto eq = equilibrium(rho_,vel); - - for(uint64_t i = 0u; i < Desc::Q; ++i){ - dfs({i}).set(dfs({i}).get() + (1.0 / relaxation_) * (eq[i] - dfs({i}).get())); - } - } -}; -} -} - -constexpr size_t dim_size = 2; -constexpr size_t dim_x = 128; -constexpr size_t dim_y = 128; - - -template -void apply_for_cells(Func&& func, saw::data& dat){ - for(std::size_t i = 0; i < dat.meta().at({1u}).get(); ++i){ - for(std::size_t j = 0; j < dat.meta().at({0u}).get(); ++j){ - saw::data di{i}; - saw::data dj{j}; - auto& cell_v = dat({{dj,di}}); - func(cell_v, j, i); - } - } -} - -void set_geometry(saw::data& latt){ - using namespace kel::lbm; - apply_for_cells([](auto& cell, std::size_t i, std::size_t j){ - uint8_t val = 0; - if(j == 1){ - val = 2u; - } - if(i == 1 || (i+2) == dim_x || (j+2) == dim_y){ - val = 3u; - } - if(i > 1 && (i+2) < dim_x && j > 1 && (j+2) < dim_y){ - val = 1u; - } - if(i == 0 || j == 0 || (i+1) == dim_x || (j+1) == dim_y){ - val = 0u; - } - cell.template get<"info">()(0u).set(val); - }, latt); -} - -void set_initial_conditions(saw::data& latt){ - using namespace kel::lbm; - - typename saw::native_data_type::type rho = 1.0; - { - std::array::type, sch::D2Q9::D> vel = {0.0, 0.0}; - auto eq = equilibrium(rho, vel); - - apply_for_cells([&eq](auto& cell, std::size_t i, std::size_t j){ - (void) i; - (void) j; - auto& dfs = cell.template get<"dfs">(); - auto& dfs_old = cell.template get<"dfs_old">(); - auto info = cell.template get<"info">()(0u).get(); - for(uint64_t k = 0; k < sch::D2Q9::Q; ++k){ - dfs(k).set(eq[k]); - dfs_old(k).set(eq[k]); - } - }, latt); - } - - { - std::array::type, sch::D2Q9::D> vel = {0.1, 0.0}; - auto eq = equilibrium(rho, vel); - - apply_for_cells([&eq](auto& cell, std::size_t i, std::size_t j){ - (void) i; - (void) j; - auto& dfs = cell.template get<"dfs">(); - auto& dfs_old = cell.template get<"dfs_old">(); - auto info = cell.template get<"info">()(0u).get(); - if(info == 2u){ - for(uint64_t k = 0; k < sch::D2Q9::Q; ++k){ - dfs(k).set(eq[k]); - dfs_old(k).set(eq[k]); - } - } - }, latt); - } -} - -void lbm_step( - saw::data& latt, - bool even_step -){ - using namespace kel::lbm; - using dfi = df_info; - - component coll; - coll.relaxation_ = 0.5384; - component rho_coll; - rho_coll.relaxation_ = 0.5384; - rho_coll.rho_ = 1.0; - - component bb; - component bb_lid; - bb_lid.lid_vel = {0.1,0.0}; - - // Collide - apply_for_cells([&](auto& cell, std::size_t i, std::size_t j){ - auto& df = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); - auto& info = cell.template get<"info">(); - - auto info_val = info({0u}).get(); - switch(info_val){ - case 1u: - coll.apply(df); - break; - case 2u: - // bb.apply(df); - bb_lid.apply(df); - break; - case 3u: - bb.apply(df); - break; - } - }, latt); - - // Stream - for(uint64_t i = 1; (i+1) < latt.template get_dim_size<0>().get(); ++i){ - for(uint64_t j = 1; (j+1) < latt.template get_dim_size<1>().get(); ++j){ - auto& cell = latt({{i,j}}); - auto& df_new = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">(); - auto& info_new = cell.template get<"info">(); - - if(info_new({0u}).get() > 0u && info_new({0u}).get() != 2u){ - for(uint64_t k = 0u; k < sch::D2Q9::Q; ++k){ - auto dir = dfi::directions[dfi::opposite_index[k]]; - auto& cell_dir_old = latt({{i+dir[0],j+dir[1]}}); - - auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">(); - auto& info_old = cell_dir_old.template get<"info">(); - - if( info_old({0}).get() == 2u ){ - auto& df_old_loc = even_step ? latt({{i,j}}).template get<"dfs_old">() : latt({{i,j}}).template get<"dfs">(); - df_new({k}) = df_old_loc({dfi::opposite_index.at(k)}) - 2.0 * dfi::inv_cs2 * dfi::weights.at(k) * 1.0 * ( bb_lid.lid_vel[0] * dir[0] + bb_lid.lid_vel[1] * dir[1]); - // dfs({dfi::opposite_index.at(i)}) = dfs_cpy({i}) - 2.0 * dfi::weights[i] * 1.0 * ( lid_vel[0] * dfi::directions[i][0] + lid_vel[1] * dfi::directions[i][1]) * dfi::inv_cs2; - } else { - df_new({k}) = df_old({k}); - } - } - } - } - } -} - -int main(){ - using namespace kel::lbm; - - saw::data> dim{{dim_x, dim_y}}; - - saw::data lattice{dim}; - - // auto& df_field = lattices.at(0).template get<"dfs">(); - //for(uint64_t i = 0; i < df_field.get_dim_size<0u>(); ++i){ - // lattices.at(i) = {dim_x, dim_y}; - //} - - /** - * Set meta information describing what this cell is - */ - set_geometry(lattice); - /** - * - */ - set_initial_conditions(lattice); - - /** - * Timeloop - */ - - /** - * Print basic setup info - */ - apply_for_cells([](auto& cell, std::size_t i, std::size_t j){ - // Not needed - (void) j; - std::cout<<(static_cast(cell.template get<"info">()({0}).get())); - if( (i+1) < dim_x){ - std::cout<<" "; - }else{ - std::cout<<"\n"; - } - }, lattice); - - uint64_t lattice_steps = 512000u; - bool even_step = true; - - uint64_t print_every = 256u; - uint64_t file_no = 0u; - - for(uint64_t step = 0; step < lattice_steps; ++step){ - - if(step % print_every == 0u){ - std::cout<<"\n"; - typename saw::native_data_type::type sum = 0.0; - - apply_for_cells([&](auto& cell, std::size_t i, std::size_t j){ - auto& dfs = cell.template get<"dfs">(); - typename saw::native_data_type::type rho; - std::array::type, sch::D2Q9::D> vel; - compute_rho_u(dfs,rho,vel); - - if(i > 1 && (i+2) < dim_x && j > 1 && (j+2) < dim_y){ - sum += rho; - } - auto angle = std::atan2(vel[1],vel[0]); - - auto dir_char = [&]() -> std::string_view { - constexpr auto pi = M_PI; - if((vel[1] * vel[1] + vel[0]*vel[0]) < 1e-16){ - return "•"; - } - if(angle > 7.0 / 8.0 * pi){ - return "←"; - } - if(angle > 5.0 / 8.0 * pi){ - return "↖"; - } - if(angle > 3.0 / 8.0 * pi){ - return "↑"; - } - if(angle > 1.0 / 8.0 * pi){ - return "↗"; - } - if(angle > -1.0 / 8.0 * pi){ - return "→"; - } - if(angle > -3.0 / 8.0 * pi){ - return "↘"; - } - if(angle > -5.0 / 8.0 * pi){ - return "↓"; - } - if(angle > -7.0 / 8.0 * pi){ - return "↙"; - } - return "←"; - }(); - std::cout<() : cell.template get<"dfs">(); - - typename saw::native_data_type::type rho; - std::array::type, sch::D2Q9::D> vel; - compute_rho_u(dfs,rho,vel); - - vtk_file << static_cast(vel[0u]) << " " << static_cast(vel[1u])<<" 0.0\n"; - }, lattice); - - - ++file_no; - } - - lbm_step(lattice, even_step); - - even_step = !even_step; - } - - /** - * Flush cout - */ - std::cout<<"\n\n"; - std::cout.flush(); - return 0; -} diff --git a/examples/SConscript b/examples/SConscript new file mode 100644 index 0000000..c9d85cb --- /dev/null +++ b/examples/SConscript @@ -0,0 +1,33 @@ +#!/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.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, ['cavity_2d.cpp'], shared=False); +examples_env.cavity_2d = examples_env.Program('#bin/cavity_2d', [env.library_static, examples_objects]); + +# Set Alias +env.examples = [examples_env.cavity_2d]; +env.Alias('examples', env.examples); + +if env["build_examples"]: + env.targets += ['examples']; + env.Install('$prefix/bin/', env.examples); +#endif diff --git a/examples/cavity_2d.cpp b/examples/cavity_2d.cpp new file mode 100644 index 0000000..3e2e064 --- /dev/null +++ b/examples/cavity_2d.cpp @@ -0,0 +1,564 @@ +#include "../c++/descriptor.hpp" + +/** + */ +#include + +#include +#include +#include + +namespace kel { +namespace lbm { +namespace sch { +using namespace saw::schema; + +/** + * Basic distribution function + * Base type + * D + * Q + * Scalar factor + * D factor + * Q factor + */ +using T = Float32; +using D2Q5 = Descriptor<2,5>; +using D2Q9 = Descriptor<2,9>; + +template +using DfCell = Cell; + +template +using CellInfo = Cell; + +/** + * Basic type for simulation + */ +template +using CellStruct = Struct< + Member, "dfs">, + Member, "dfs_old">, + Member, "info"> +>; + + +using CavityFieldD2Q9 = CellField>; +} + +/** + * Calculate the macroscopic variables rho and u in Lattice Units. + */ +template +void compute_rho_u ( + saw::data>& dfs, + typename saw::native_data_type::type& rho, + std::array::type, 2>& vel + ) +{ + using dfi = df_info; + + rho = 0; + std::fill(vel.begin(), vel.end(), 0); + + for(size_t i = 0; i < Desc::Q; ++i){ + rho += dfs(i).get(); + vel[0] += dfi::directions[i][0] * dfs(i).get(); + vel[1] += dfi::directions[i][1] * dfs(i).get(); + } + + vel[0] /= rho; + vel[1] /= rho; +} + +/** + * Unsure if feasible. I would have a rho normalization on the dfs and then I would use the const rho_computation + */ +template +void compute_const_rho_u ( + saw::data>& dfs, + const typename saw::native_data_type::type rho, + std::array::type, 2>& vel + ) +{ + using dfi = df_info; + + saw::native_data_type::type real_rho = 0; + std::fill(vel.begin(), vel.end(), 0); + + for(size_t i = 0; i < Desc::Q; ++i){ + real_rho += dfs(i).get(); + vel[0] += dfi::directions[i][0] * dfs(i).get(); + vel[1] += dfi::directions[i][1] * dfs(i).get(); + } + for(size_t i = 0; i < Desc::Q; ++i){ + dfs(i).set(dfs(i).get() * (rho/real_rho)); + } + + vel[0] *= real_rho / (rho*rho); + vel[1] *= real_rho / (rho*rho); +} + +/** + * Calculates the equilibrium for each direction + */ +template +std::array::type,Desc::Q> equilibrium( + typename saw::native_data_type::type rho, + const std::array::type, Desc::D>& vel +){ + using dfi = df_info; + + typename std::array::type,Desc::Q> eq; + + for(std::size_t i = 0; i < eq.size(); ++i){ + auto vel_c = (vel[0]*dfi::directions[i][0] + vel[1]*dfi::directions[i][1]); + auto vel_c_cs2 = vel_c * dfi::inv_cs2; + eq[i] = dfi::weights[i] * rho * ( + 1.0 + + vel_c_cs2 + - dfi::inv_cs2 * 0.5 * ( vel[0] * vel[0] + vel[1] * vel[1] ) + + vel_c_cs2 * vel_c_cs2 * 0.5 + ); + } + + return eq; +} + +/** + * A reason for why a component based design is really good can be seen in my LR solver example + * + * Add Expression Templates and you're golden. + */ +template +class component; + +/* +template +class df_cell_view; +*/ +/** + * Minor helper for the AA-Pull Pattern, so I can use only one lattice + * + * Am I sure I want to use AA this way? + * Esoteric Twist technically reduces the needed memory access footprint + */ +/* +template +class df_cell_view, Encode> { +public: + using Schema = sch::Cell; +private: + std::array::type>*, QN> view_; +public: + df_cell_view(const std::array::type>*, QN>& view): + view_{view} + {} +}; +*/ +namespace cmpt { +struct BounceBack{}; +struct MovingWall {}; +struct BGK {}; +struct ConstRhoBGK {}; +} + + +/** + * Full-Way BounceBack. I suspect that the moving wall requires half-way bounce back. + */ +template +class component { +public: + + void apply(saw::data>& dfs){ + using dfi = df_info; + + // Technically use .copy() + auto df_cpy = dfs; + + for(uint64_t i = 1u; i < Desc::Q; ++i){ + dfs({i}) = df_cpy({dfi::opposite_index.at(i)}); + } + } +}; + +/** + * Full-Way moving wall Bounce back, something is not right here. + * Technically it should reflect properly. + */ +template +class component { +public: + std::array::type, Desc::D> lid_vel; + +public: + void apply( + saw::data>& dfs + ){ + using dfi = df_info; + + // Technically use .copy() + /* + auto dfs_cpy = dfs; + + for(uint64_t i = 0u; i < Desc::Q; ++i){ + dfs({dfi::opposite_index.at(i)}) = dfs_cpy({i}) - 2.0 * dfi::weights[i] * 1.0 * ( lid_vel[0] * dfi::directions[i][0] + lid_vel[1] * dfi::directions[i][1]) * dfi::inv_cs2; + } + */ + } +}; + +template +class component { +public: + typename saw::native_data_type::type relaxation_; +public: + void apply(saw::data>& dfs){ + typename saw::native_data_type::type rho; + std::array::type, Desc::D> vel; + compute_rho_u(dfs,rho,vel); + auto eq = equilibrium(rho,vel); + + for(uint64_t i = 0u; i < Desc::Q; ++i){ + dfs({i}).set(dfs({i}).get() + (1.0 / relaxation_) * (eq[i] - dfs({i}).get())); + } + } +}; + +template +class component { +public: + typename saw::native_data_type::type relaxation_; + typename saw::native_data_type::type rho_; +public: + void apply(saw::data>& dfs){ + std::array::type, Desc::D> vel; + compute_const_rho_u(dfs,rho_,vel); + auto eq = equilibrium(rho_,vel); + + for(uint64_t i = 0u; i < Desc::Q; ++i){ + dfs({i}).set(dfs({i}).get() + (1.0 / relaxation_) * (eq[i] - dfs({i}).get())); + } + } +}; +} +} + +constexpr size_t dim_size = 2; +constexpr size_t dim_x = 128; +constexpr size_t dim_y = 128; + + +template +void apply_for_cells(Func&& func, saw::data& dat){ + for(std::size_t i = 0; i < dat.meta().at({1u}).get(); ++i){ + for(std::size_t j = 0; j < dat.meta().at({0u}).get(); ++j){ + saw::data di{i}; + saw::data dj{j}; + auto& cell_v = dat({{dj,di}}); + func(cell_v, j, i); + } + } +} + +void set_geometry(saw::data& latt){ + using namespace kel::lbm; + apply_for_cells([](auto& cell, std::size_t i, std::size_t j){ + uint8_t val = 0; + if(j == 1){ + val = 2u; + } + if(i == 1 || (i+2) == dim_x || (j+2) == dim_y){ + val = 3u; + } + if(i > 1 && (i+2) < dim_x && j > 1 && (j+2) < dim_y){ + val = 1u; + } + if(i == 0 || j == 0 || (i+1) == dim_x || (j+1) == dim_y){ + val = 0u; + } + cell.template get<"info">()(0u).set(val); + }, latt); +} + +void set_initial_conditions(saw::data& latt){ + using namespace kel::lbm; + + typename saw::native_data_type::type rho = 1.0; + { + std::array::type, sch::D2Q9::D> vel = {0.0, 0.0}; + auto eq = equilibrium(rho, vel); + + apply_for_cells([&eq](auto& cell, std::size_t i, std::size_t j){ + (void) i; + (void) j; + auto& dfs = cell.template get<"dfs">(); + auto& dfs_old = cell.template get<"dfs_old">(); + auto info = cell.template get<"info">()(0u).get(); + for(uint64_t k = 0; k < sch::D2Q9::Q; ++k){ + dfs(k).set(eq[k]); + dfs_old(k).set(eq[k]); + } + }, latt); + } + + { + std::array::type, sch::D2Q9::D> vel = {0.1, 0.0}; + auto eq = equilibrium(rho, vel); + + apply_for_cells([&eq](auto& cell, std::size_t i, std::size_t j){ + (void) i; + (void) j; + auto& dfs = cell.template get<"dfs">(); + auto& dfs_old = cell.template get<"dfs_old">(); + auto info = cell.template get<"info">()(0u).get(); + if(info == 2u){ + for(uint64_t k = 0; k < sch::D2Q9::Q; ++k){ + dfs(k).set(eq[k]); + dfs_old(k).set(eq[k]); + } + } + }, latt); + } +} + +void lbm_step( + saw::data& latt, + bool even_step +){ + using namespace kel::lbm; + using dfi = df_info; + + component coll; + coll.relaxation_ = 0.5384; + component rho_coll; + rho_coll.relaxation_ = 0.5384; + rho_coll.rho_ = 1.0; + + component bb; + component bb_lid; + bb_lid.lid_vel = {0.1,0.0}; + + // Collide + apply_for_cells([&](auto& cell, std::size_t i, std::size_t j){ + auto& df = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + auto& info = cell.template get<"info">(); + + auto info_val = info({0u}).get(); + switch(info_val){ + case 1u: + coll.apply(df); + break; + case 2u: + // bb.apply(df); + bb_lid.apply(df); + break; + case 3u: + bb.apply(df); + break; + } + }, latt); + + // Stream + for(uint64_t i = 1; (i+1) < latt.template get_dim_size<0>().get(); ++i){ + for(uint64_t j = 1; (j+1) < latt.template get_dim_size<1>().get(); ++j){ + auto& cell = latt({{i,j}}); + auto& df_new = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">(); + auto& info_new = cell.template get<"info">(); + + if(info_new({0u}).get() > 0u && info_new({0u}).get() != 2u){ + for(uint64_t k = 0u; k < sch::D2Q9::Q; ++k){ + auto dir = dfi::directions[dfi::opposite_index[k]]; + auto& cell_dir_old = latt({{i+dir[0],j+dir[1]}}); + + auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">(); + auto& info_old = cell_dir_old.template get<"info">(); + + if( info_old({0}).get() == 2u ){ + auto& df_old_loc = even_step ? latt({{i,j}}).template get<"dfs_old">() : latt({{i,j}}).template get<"dfs">(); + df_new({k}) = df_old_loc({dfi::opposite_index.at(k)}) - 2.0 * dfi::inv_cs2 * dfi::weights.at(k) * 1.0 * ( bb_lid.lid_vel[0] * dir[0] + bb_lid.lid_vel[1] * dir[1]); + // dfs({dfi::opposite_index.at(i)}) = dfs_cpy({i}) - 2.0 * dfi::weights[i] * 1.0 * ( lid_vel[0] * dfi::directions[i][0] + lid_vel[1] * dfi::directions[i][1]) * dfi::inv_cs2; + } else { + df_new({k}) = df_old({k}); + } + } + } + } + } +} + +int main(){ + using namespace kel::lbm; + + saw::data> dim{{dim_x, dim_y}}; + + saw::data lattice{dim}; + + // auto& df_field = lattices.at(0).template get<"dfs">(); + //for(uint64_t i = 0; i < df_field.get_dim_size<0u>(); ++i){ + // lattices.at(i) = {dim_x, dim_y}; + //} + + /** + * Set meta information describing what this cell is + */ + set_geometry(lattice); + /** + * + */ + set_initial_conditions(lattice); + + /** + * Timeloop + */ + + /** + * Print basic setup info + */ + apply_for_cells([](auto& cell, std::size_t i, std::size_t j){ + // Not needed + (void) j; + std::cout<<(static_cast(cell.template get<"info">()({0}).get())); + if( (i+1) < dim_x){ + std::cout<<" "; + }else{ + std::cout<<"\n"; + } + }, lattice); + + uint64_t lattice_steps = 512000u; + bool even_step = true; + + uint64_t print_every = 256u; + uint64_t file_no = 0u; + + for(uint64_t step = 0; step < lattice_steps; ++step){ + + if(step % print_every == 0u){ + std::cout<<"\n"; + typename saw::native_data_type::type sum = 0.0; + + + apply_for_cells([&](auto& cell, std::size_t i, std::size_t j){ + auto& dfs = cell.template get<"dfs">(); + typename saw::native_data_type::type rho; + std::array::type, sch::D2Q9::D> vel; + compute_rho_u(dfs,rho,vel); + + if(i > 1 && (i+2) < dim_x && j > 1 && (j+2) < dim_y){ + sum += rho; + } + auto angle = std::atan2(vel[1],vel[0]); + double vel_mag = vel[1] * vel[1] + vel[0] * vel[0]; + + auto dir_char = [&]() -> std::string_view { + /** + * Flipped y due to print order + */ + constexpr auto pi = M_PI; + return "■"; + if(vel_mag < 1e-4){ + return "•"; + } + if(angle > 7.0 / 8.0 * pi){ + return "←"; + } + if(angle > 5.0 / 8.0 * pi){ + return "↙"; + } + if(angle > 3.0 / 8.0 * pi){ + return "↓"; + } + if(angle > 1.0 / 8.0 * pi){ + return "↘"; + } + if(angle > -1.0 / 8.0 * pi){ + return "→"; + } + if(angle > -3.0 / 8.0 * pi){ + return "↗"; + } + if(angle > -5.0 / 8.0 * pi){ + return "↑"; + } + if(angle > -7.0 / 8.0 * pi){ + return "↖"; + } + return "←"; + }(); + + std::array rgb_start{64,64,255}; + std::array rgb_stop{255,64,64}; + std::array rgb_middle{255,255,255}; + + + auto col_interpol = [&](){ + std::array rgb_interpol = rgb_start; + double vel_mag_cut = std::min(1.0,std::max(vel_mag/(0.07*0.07),0.0)); + + if(vel_mag_cut < 0.5){ + uint32_t vel_mag_i = static_cast(2.0 * vel_mag_cut * 256); + for(uint8_t i = 0u; i < 3u; ++i){ + rgb_interpol[i] = (rgb_middle[i] * vel_mag_i + (256-vel_mag_i) * rgb_start[i]) / 256; + } + }else{ + uint32_t vel_mag_i = static_cast((2.0*(vel_mag_cut-0.5)) * 256); + for(uint8_t i = 0u; i < 3u; ++i){ + rgb_interpol[i] = (rgb_stop[i] * vel_mag_i + (256-vel_mag_i) * rgb_middle[i]) / 256; + } + } + return rgb_interpol; + }; + auto rgb_interpol = col_interpol(); + + std::cout<<"\x1b[38;2;"<() : cell.template get<"dfs">(); + + typename saw::native_data_type::type rho; + std::array::type, sch::D2Q9::D> vel; + compute_rho_u(dfs,rho,vel); + + vtk_file << static_cast(vel[0u]) << " " << static_cast(vel[1u])<<" 0.0\n"; + }, lattice); + + + ++file_no; + } + + lbm_step(lattice, even_step); + + even_step = !even_step; + } + + /** + * Flush cout + */ + std::cout<<"\n\n"; + std::cout.flush(); + return 0; +} -- cgit v1.2.3